---
title: "Replication Stage 2 Report"
author: "Adeline Lo"
date: "March 2021"
output: 
  html_document: 
    code_folding: hide
    pandoc_args: ["--lua-filter=color-text.lua"]
    toc: true
    toc_float: true
    number_sections: true
    theme: paper
  pdf_document: 
    pandoc_args: ["--lua-filter=color-text.lua"]
    keep_tex: true
---
<style type="text/css">
.main-container {
  max-width: 1800px;
  margin-left: auto;
  margin-right: auto;
}
</style>

```{cat, engine.opts = list(file = "color-text.lua")}
Span = function(span)
  color = span.attributes['color']
  -- if no color attribute, return unchange
  if color == nil then return span end
  
  -- tranform to <span style="color: red;"></span>
  if FORMAT:match 'html' then
    -- remove color attributes
    span.attributes['color'] = nil
    -- use style attribute instead
    span.attributes['style'] = 'color: ' .. color .. ';'
    -- return full span element
    return span
  elseif FORMAT:match 'latex' then
    -- remove color attributes
    span.attributes['color'] = nil
    -- encapsulate in latex code
    table.insert(
      span.content, 1,
      pandoc.RawInline('latex', '\\textcolor{'..color..'}{')
    )
    table.insert(
      span.content,
      pandoc.RawInline('latex', '}')
    )
    -- returns only span content
    return span.content
  else
    -- for other format return unchanged
    return span
  end
end
```

```{css, echo=FALSE}
.pretty {
  background-color: white;
  border: 3px solid blue;
  font-weight: bold;
  font-family: Courier;
}
caption {
      color: darkblue;
      font-weight: bold;
      font-size: 1.0em;
    } 
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(anytime)
library(dplyr)
library(knitr)
library(tidyverse)
#graphics
library(viridis)
library(viridisLite)
library(ggplot2)
library(tidyverse)
library(hrbrthemes)
library(ggsci)
library(xtable)
library(gridExtra)
#testing
library(sandwich)
library(lmtest)
library(rstanarm)
library(bayestestR)
library(logspline)
library(effectsize)
library(see)
#html tables
library(kableExtra)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(texreg)
#output tables
library(stargazer)
#other
library(magrittr)
library(stringr)
#design simulations
library(DeclareDesign)
library(DesignLibrary)
```


# Data cleaning

Data pulled from `Data/Mustafa-Experiment-Oct-Final.png`, `Data/Kelli-Experiment-Oct-Final.png`. Combined with`Data/Facebook Likes and Comments`, outputs into `Data` folder:

* `mustafa.rds`
* `kelli.rds`
* `df.rds`
* `test_df`: includes t-test hypotheses info of all ads

```{r Data cleaning,eval=TRUE,message=FALSE,warning=FALSE}
mustafa<-data.frame(
  treatment_name=c("COVID-US","COVID-PA","COVID-Lancaster","COVID-No-Place","No-COVID-US")#treatment name
  ,treatment=c(1,2,3,4,5) #treatment
  ,results=c(471,542,561,484,489) #link clicks
  ,reach=c(24536,27520,27191,25112,26008)
  ,impressions=c(48934,50107,50583,48987,50098)
  ,n_reactions=c(113,236,242,101,135) #combined emojis
  ,n_like=c(101,191,209,91,121)
  ,n_love=c(11,43,30,8,12)
  ,n_care=c(0,2,0,2,0)
  ,n_haha=c(0,0,3,0,0)
  ,n_wow=c(1,0,0,0,0)
  ,n_anger=c(0,0,0,0,2)
  ,n_sad=rep(0,5)
  ,n_shares=c(7,25,19,9,7)
  ,n_comments=c(11,7,9,3,7))
mustafa$prop_click<-mustafa$results/mustafa$reach
mustafa$prop_click2<-mustafa$results/mustafa$impressions

#Table 1: Mustafa part
kable(mustafa,caption="Mustafa Ad data")%>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"),full_width=F)
saveRDS(mustafa,file="Data/mustafa.rds")

mustafa_comments<-data.frame(treatment_name=rep(c("COVID-US","COVID-PA","COVID-Lancaster","COVID-No-Place","No-COVID-US")
                                                  ,times=mustafa$n_comments)#treatment name
  ,comments=c(c("Amen","Thank you","Thank you very much!","Thank you for your help.","Thank you","thx to all"
                ,"Thank you to all & God Bless all","Or bring it to us","Thank you young man!! Virtual hug"
                ,"Thank you so much","Ameen") #covid-us T1; no emoticons/images coded, only text (or text in image).
             ,c("Thank you!","You are a great man, Mustafa!","Awesome!","God Bless"
                ,"Now we can call them refugees O boy I call them illegals the people can't afford to buy food there's no work who's paying for them yes we are","Refugees are usually brought here by religious groups. They are not illegal!!!","")#covid-pa T2;
             ,c("Thank you from a person (if not handicapped) would gladly help.","Great","Pretty tricky camera work"
                ,"God Bless","I will","I believe this man is the same you know it doesn't matter what skin color we are we all believe the same it's the government having us fighting against one each other we should all get together and fight the government","Coronahoax"
                ,"Be better if they were to go back and do it in their own country"
                ,"Why aren't any refugees in non white countries? Why are they always in white countries and then the filth bags have the nerve to call us racist")#covid-la T3
             ,c("Thank you for showing kindness.","Nice, now pick up the milk!","Is it just me or does he look like  a tiny guy standing in a cart?")#covid-noplace T4
             ,c("Keep up the good work","How about just, people helping people","At first, I thought this picture was a tiny man standing in a cart","Nothing is dangerous than incomplete knowledge","how naive of you guys!!"
                ,"Deport all refugees unless thet are legal","Thanks for the show.")#no-covid-us T5
  )#end comments
)
#code tone of comments, -1 negative, 0 neutral, 1 positive
mustafa_comments$comments_tone<-c(1,1,1,1,1,1,1,-1,1,1
                                  ,1,1,1,1,1,-1,1,0,1,1
                                  ,-1,1,1,1,-1,-1,-1,1,-1,0
                                  ,1,0,0,0,-1,-1,-1)

kelli<-data.frame(
  treatment_name=c("COVID-Refugee","COVID-Immigrant","COVID-Neither","No-COVID-Refugee","No-COVID-Neither")#treatment name
  ,treatment=c(1,2,3,4,5) #treatment
  ,results=c(1261,1343,1193,1295,1121) #link clicks
  ,reach=c(44224,48320,46679,44384,46088)#number treated
  ,impressions=c(58632,63272,59221,59380,58553)#times
  ,n_reactions=c(291,534,758,208,347) #combined emojis
  ,n_like=c(278,477,678,192,330)
  ,n_love=c(8,43,67,13,14)
  ,n_care=c(0,14,12,3,3)
  ,n_haha=rep(0,5)
  ,n_wow=rep(0,5)
  ,n_anger=rep(0,5)
  ,n_sad=c(5,0,0,0,0)
  ,n_shares=c(29,43,56,16,31)
  ,n_comments=c(18,35,32,8,18))
kelli$prop_click<-kelli$results/kelli$reach
kelli$prop_click2<-kelli$results/kelli$impressions

#Table 1: Kelli part
kable(kelli,caption="Kelli Ad data")%>%
  kable_styling(bootstrap_options=c("striped","hover","condensed"),full_width=F)
saveRDS(kelli,file="Data/kelli.rds")

mdf<-data.frame(
  ad="Mustafa"
  ,treatment_name=rep(mustafa$treatment_name,times=mustafa$reach)
  ,treatment=rep(mustafa$treatment,times=mustafa$reach)
  ,click=c(rep(1,mustafa$results[1]),rep(0,(mustafa$reach[1]-mustafa$results[1]))#t1
          ,rep(1,mustafa$results[2]),rep(0,(mustafa$reach[2]-mustafa$results[2]))#t2
          ,rep(1,mustafa$results[3]),rep(0,(mustafa$reach[3]-mustafa$results[3]))#t3
          ,rep(1,mustafa$results[4]),rep(0,(mustafa$reach[4]-mustafa$results[4]))#t4
          ,rep(1,mustafa$results[5]),rep(0,(mustafa$reach[5]-mustafa$results[5])))#t5
  ,reaction=c(rep(1,mustafa$n_reactions[1]),rep(0,(mustafa$reach[1]-mustafa$n_reactions[1]))#t1
          ,rep(1,mustafa$n_reactions[2]),rep(0,(mustafa$reach[2]-mustafa$n_reactions[2]))#t2
          ,rep(1,mustafa$n_reactions[3]),rep(0,(mustafa$reach[3]-mustafa$n_reactions[3]))#t3
          ,rep(1,mustafa$n_reactions[4]),rep(0,(mustafa$reach[4]-mustafa$n_reactions[4]))#t4
          ,rep(1,mustafa$n_reactions[5]),rep(0,(mustafa$reach[5]-mustafa$n_reactions[5])))#t5

  ,like=c(rep(1,mustafa$n_like[1]),rep(0,(mustafa$reach[1]-mustafa$n_like[1]))#t1
          ,rep(1,mustafa$n_like[2]),rep(0,(mustafa$reach[2]-mustafa$n_like[2]))#t2
          ,rep(1,mustafa$n_like[3]),rep(0,(mustafa$reach[3]-mustafa$n_like[3]))#t3
          ,rep(1,mustafa$n_like[4]),rep(0,(mustafa$reach[4]-mustafa$n_like[4]))#t4
          ,rep(1,mustafa$n_like[5]),rep(0,(mustafa$reach[5]-mustafa$n_like[5])))#t5

  ,love=c(rep(1,mustafa$n_love[1]),rep(0,(mustafa$reach[1]-mustafa$n_love[1]))#t1
          ,rep(1,mustafa$n_love[2]),rep(0,(mustafa$reach[2]-mustafa$n_love[2]))#t2
          ,rep(1,mustafa$n_love[3]),rep(0,(mustafa$reach[3]-mustafa$n_love[3]))#t3
          ,rep(1,mustafa$n_love[4]),rep(0,(mustafa$reach[4]-mustafa$n_love[4]))#t4
          ,rep(1,mustafa$n_love[5]),rep(0,(mustafa$reach[5]-mustafa$n_love[5])))#t5
  ,care=c(rep(1,mustafa$n_care[1]),rep(0,(mustafa$reach[1]-mustafa$n_care[1]))#t1
          ,rep(1,mustafa$n_care[2]),rep(0,(mustafa$reach[2]-mustafa$n_care[2]))#t2
          ,rep(1,mustafa$n_care[3]),rep(0,(mustafa$reach[3]-mustafa$n_care[3]))#t3
          ,rep(1,mustafa$n_care[4]),rep(0,(mustafa$reach[4]-mustafa$n_care[4]))#t4
          ,rep(1,mustafa$n_care[5]),rep(0,(mustafa$reach[5]-mustafa$n_care[5])))#t5

  ,positive=c(rep(1,c(mustafa$n_like[1]+mustafa$n_love[1]+mustafa$n_care[1])),rep(0,(mustafa$reach[1]-c(mustafa$n_like[1]+mustafa$n_love[1]+mustafa$n_care[1])))#t1
    ,rep(1,c(mustafa$n_like[2]+mustafa$n_love[2]+mustafa$n_care[2])),rep(0,(mustafa$reach[2]-c(mustafa$n_like[2]+mustafa$n_love[2]+mustafa$n_care[2])))#t2
          ,rep(1,c(mustafa$n_like[3]+mustafa$n_love[3]+mustafa$n_care[3])),rep(0,(mustafa$reach[3]-c(mustafa$n_like[3]+mustafa$n_love[3]+mustafa$n_care[3])))#t3
          ,rep(1,c(mustafa$n_like[4]+mustafa$n_love[4]+mustafa$n_care[4])),rep(0,(mustafa$reach[4]-c(mustafa$n_like[4]+mustafa$n_love[4]+mustafa$n_care[4])))#t4
          ,rep(1,c(mustafa$n_like[5]+mustafa$n_love[5]+mustafa$n_care[5])),rep(0,(mustafa$reach[5]-c(mustafa$n_like[5]+mustafa$n_love[5]+mustafa$n_care[5]))))#t5
,negative=c(rep(1,c(mustafa$n_anger[1]+mustafa$n_sad[1])),rep(0,(mustafa$reach[1]-c(mustafa$n_anger[1]+mustafa$n_sad[1])))#t1
          ,rep(1,c(mustafa$n_anger[2]+mustafa$n_sad[2])),rep(0,(mustafa$reach[2]-c(mustafa$n_anger[2]+mustafa$n_sad[2])))#t2
          ,rep(1,c(mustafa$n_anger[3]+mustafa$n_sad[3])),rep(0,(mustafa$reach[3]-c(mustafa$n_anger[3]+mustafa$n_sad[3])))#t3
          ,rep(1,c(mustafa$n_anger[4]+mustafa$n_sad[4])),rep(0,(mustafa$reach[4]-c(mustafa$n_anger[4]+mustafa$n_sad[4])))#t4
          ,rep(1,c(mustafa$n_anger[5]+mustafa$n_sad[5])),rep(0,(mustafa$reach[5]-c(mustafa$n_anger[5]+mustafa$n_sad[5]))))#t5
)

kdf<-data.frame(
  ad="Kelli"
  ,treatment_name=rep(kelli$treatment_name,times=kelli$reach)
  ,treatment=rep(kelli$treatment,times=kelli$reach)
  ,click=c(rep(1,kelli$results[1]),rep(0,(kelli$reach[1]-kelli$results[1]))#t1
          ,rep(1,kelli$results[2]),rep(0,(kelli$reach[2]-kelli$results[2]))#t2
          ,rep(1,kelli$results[3]),rep(0,(kelli$reach[3]-kelli$results[3]))#t3
          ,rep(1,kelli$results[4]),rep(0,(kelli$reach[4]-kelli$results[4]))#t4
          ,rep(1,kelli$results[5]),rep(0,(kelli$reach[5]-kelli$results[5]))#t5
          )
  ,reaction=c(rep(1,kelli$n_reactions[1]),rep(0,(kelli$reach[1]-kelli$n_reactions[1]))#t1
          ,rep(1,kelli$n_reactions[2]),rep(0,(kelli$reach[2]-kelli$n_reactions[2]))#t2
          ,rep(1,kelli$n_reactions[3]),rep(0,(kelli$reach[3]-kelli$n_reactions[3]))#t3
          ,rep(1,kelli$n_reactions[4]),rep(0,(kelli$reach[4]-kelli$n_reactions[4]))#t4
          ,rep(1,kelli$n_reactions[5]),rep(0,(kelli$reach[5]-kelli$n_reactions[5])))#t5

  ,like=c(rep(1,kelli$n_like[1]),rep(0,(kelli$reach[1]-kelli$n_like[1]))#t1
          ,rep(1,kelli$n_like[2]),rep(0,(kelli$reach[2]-kelli$n_like[2]))#t2
          ,rep(1,kelli$n_like[3]),rep(0,(kelli$reach[3]-kelli$n_like[3]))#t3
          ,rep(1,kelli$n_like[4]),rep(0,(kelli$reach[4]-kelli$n_like[4]))#t4
          ,rep(1,kelli$n_like[5]),rep(0,(kelli$reach[5]-kelli$n_like[5])))#t5

  ,love=c(rep(1,kelli$n_love[1]),rep(0,(kelli$reach[1]-kelli$n_love[1]))#t1
          ,rep(1,kelli$n_love[2]),rep(0,(kelli$reach[2]-kelli$n_love[2]))#t2
          ,rep(1,kelli$n_love[3]),rep(0,(kelli$reach[3]-kelli$n_love[3]))#t3
          ,rep(1,kelli$n_love[4]),rep(0,(kelli$reach[4]-kelli$n_love[4]))#t4
          ,rep(1,kelli$n_love[5]),rep(0,(kelli$reach[5]-kelli$n_love[5])))#t5
  ,care=c(rep(1,kelli$n_care[1]),rep(0,(kelli$reach[1]-kelli$n_care[1]))#t1
          ,rep(1,kelli$n_care[2]),rep(0,(kelli$reach[2]-kelli$n_care[2]))#t2
          ,rep(1,kelli$n_care[3]),rep(0,(kelli$reach[3]-kelli$n_care[3]))#t3
          ,rep(1,kelli$n_care[4]),rep(0,(kelli$reach[4]-kelli$n_care[4]))#t4
          ,rep(1,kelli$n_care[5]),rep(0,(kelli$reach[5]-kelli$n_care[5])))#t5

  ,positive=c(rep(1,c(kelli$n_like[1]+kelli$n_love[1]+kelli$n_care[1])),rep(0,(kelli$reach[1]-c(kelli$n_like[1]+kelli$n_love[1]+kelli$n_care[1])))#t1
    ,rep(1,c(kelli$n_like[2]+kelli$n_love[2]+kelli$n_care[2])),rep(0,(kelli$reach[2]-c(kelli$n_like[2]+kelli$n_love[2]+kelli$n_care[2])))#t2
          ,rep(1,c(kelli$n_like[3]+kelli$n_love[3]+kelli$n_care[3])),rep(0,(kelli$reach[3]-c(kelli$n_like[3]+kelli$n_love[3]+kelli$n_care[3])))#t3
          ,rep(1,c(kelli$n_like[4]+kelli$n_love[4]+kelli$n_care[4])),rep(0,(kelli$reach[4]-c(kelli$n_like[4]+kelli$n_love[4]+kelli$n_care[4])))#t4
          ,rep(1,c(kelli$n_like[5]+kelli$n_love[5]+kelli$n_care[5])),rep(0,(kelli$reach[5]-c(kelli$n_like[5]+kelli$n_love[5]+kelli$n_care[5]))))#t5
,negative=c(rep(1,c(kelli$n_anger[1]+kelli$n_sad[1])),rep(0,(kelli$reach[1]-c(kelli$n_anger[1]+kelli$n_sad[1])))#t1
          ,rep(1,c(kelli$n_anger[2]+kelli$n_sad[2])),rep(0,(kelli$reach[2]-c(kelli$n_anger[2]+kelli$n_sad[2])))#t2
          ,rep(1,c(kelli$n_anger[3]+kelli$n_sad[3])),rep(0,(kelli$reach[3]-c(kelli$n_anger[3]+kelli$n_sad[3])))#t3
          ,rep(1,c(kelli$n_anger[4]+kelli$n_sad[4])),rep(0,(kelli$reach[4]-c(kelli$n_anger[4]+kelli$n_sad[4])))#t4
          ,rep(1,c(kelli$n_anger[5]+kelli$n_sad[5])),rep(0,(kelli$reach[5]-c(kelli$n_anger[5]+kelli$n_sad[5]))))#t5
)

df <- rbind(mdf,kdf)

test_df<-data.frame(
  Hypothesis=c(rep("H1 Covid effect",3), rep("H2 Location effect",6),rep("H3 Refugee effect",3))
  ,Group1=NA
  ,Group2=NA
  ,Group1Mean=NA
  ,Group2Mean=NA
  ,Mean=NA
  ,Upper=NA
  ,Lower=NA
  ,p=NA
)#holds all ttest outputs
```

# Main Hypotheses

## H1: Covid or No Covid
```{r TestingH1,eval=TRUE,warning=FALSE,message=FALSE}
### Covid US vs No Covid US [T1-T5]
#**Difference between Mustafa T1 and Mustafa T5**
mdf15<-subset(mdf,treatment==1|treatment==5)
mdf15$treatment_name<-as.factor(mdf15$treatment_name)
mdf15$treatment_name<-relevel(mdf15$treatment_name, ref="COVID-US")
tmp<-t.test(click~treatment_name,mdf15)
test_df[1,c("Group1","Group2")]<-levels(mdf15$treatment_name)
test_df[1,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[1]<-tmp$estimate[1]-tmp$estimate[2]
test_df[1,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[1]<-"Mustafa"
test_df$p[1]<-tmp$p.value
### Covid Refugee vs No-Covid Refugee [T1-T4]
#**Difference between Kelli T1 and Kelli T4**
kdf14<-subset(kdf,treatment==1|treatment==4)
kdf14$treatment_name<-as.factor(kdf14$treatment_name)
kdf14$treatment_name<-relevel(kdf14$treatment_name, ref="COVID-Refugee")
tmp<-t.test(click~treatment_name,kdf14)
test_df[2,c("Group1","Group2")]<-levels(kdf14$treatment_name)
test_df[2,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[2]<-tmp$estimate[1]-tmp$estimate[2]
test_df[2,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[2]<-"Kelli"
test_df$p[2]<-tmp$p.value
### Covid Neither vs No-Covid Neither [T3-T5]
#**Difference between Kelli T3 and Kelli T5**
kdf35<-subset(kdf,treatment==3|treatment==5)
kdf35$treatment_name<-as.factor(kdf35$treatment_name)
kdf35$treatment_name<-relevel(kdf35$treatment_name, ref="COVID-Neither")
tmp<-t.test(click~treatment_name,kdf35)
test_df[3,c("Group1","Group2")]<-levels(kdf35$treatment_name)
test_df[3,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[3]<-tmp$estimate[1]-tmp$estimate[2]
test_df[3,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[3]<-"Kelli"
test_df$p[3]<-tmp$p.value
```

```{r,eval=FALSE}
### Covid US vs No Covid US [T1-T5]
#**Difference between Mustafa T1 and Mustafa T5**
#**Bayes' factor**
#Given the observed experimental data, has the null hypothesis of no difference between No-Covid-US and Covid-US or absence of an effect become more, or less credible?
#We assign informative priors assuming a linear regression model:
#$click_i \sim \text{Normal} (\alpha + \beta Treatment_i, \sigma)$
#where we take into account of the low click rates for FB ads generally as well as our assumed treatment effect sizes (small around 0.1) and consider that approximately $\beta\in(-0.5,0.5)$. We therefore set the prior mean at the midpoint of the interval and allow for some variation around either side ($\sigma=0.2$).
covid_prior_small<- normal(location = 0, scale = 0.02)
covid_prior_med <- normal(location = 0, scale = 0.04)
covid_prior_large <- normal(location = 0, scale = 0.16)
mod_small <- stan_glm(click ~ treatment_name, data = mdf15, prior = covid_prior_small)
mod_med <- stan_glm(click ~ treatment_name, data = mdf15, prior = covid_prior_med)
mod_large <- stan_glm(click ~ treatment_name, data = mdf15, prior = covid_prior_large)
##compute bayes factor as change from prior odds to posterior odds, then save:
mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
saveRDS(mod_bf_small,file="Data/mod_bf_small_M1-M5.rds")
mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
saveRDS(mod_bf_med,file="Data/mod_bf_med_M1-M5.rds")
mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
saveRDS(mod_bf_large,file="Data/mod_bf_large_M1-M5.rds")
### Covid Refugee vs No-Covid Refugee [T1-T4]
#**Difference between Kelli T1 and Kelli T4**
#**Bayes' factor**
covid_prior_small<- normal(location = 0, scale = 0.02)
covid_prior_med <- normal(location = 0, scale = 0.04)
covid_prior_large <- normal(location = 0, scale = 0.16)
mod_small <- stan_glm(click ~ treatment_name, data = kdf14, prior = covid_prior_small)
mod_med <- stan_glm(click ~ treatment_name, data = kdf14, prior = covid_prior_med)
mod_large <- stan_glm(click ~ treatment_name, data = kdf14, prior = covid_prior_large)
##compute bayes factor as change from prior odds to posterior odds:
mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
saveRDS(mod_bf_small,file="Data/mod_bf_small_K1-K4.rds")
mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
saveRDS(mod_bf_med,file="Data/mod_bf_med_K1-K4.rds")
mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
saveRDS(mod_bf_large,file="Data/mod_bf_large_K1-K4.rds")

### Covid Neither vs No-Covid Neither [T3-T5]
#**Difference between Kelli T3 and Kelli T5**
#**Bayes' factor**
covid_prior_small<- normal(location = 0, scale = 0.02)
covid_prior_med <- normal(location = 0, scale = 0.04)
covid_prior_large <- normal(location = 0, scale = 0.16)
mod_small <- stan_glm(click ~ treatment_name, data = kdf35, prior = covid_prior_small)
mod_med <- stan_glm(click ~ treatment_name, data = kdf35, prior = covid_prior_med)
mod_large <- stan_glm(click ~ treatment_name, data = kdf35, prior = covid_prior_large)
##compute bayes factor as change from prior odds to posterior odds:
mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
saveRDS(mod_bf_small,file="Data/mod_bf_small_K3-K5.rds") 
mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
saveRDS(mod_bf_med,file="Data/mod_bf_med_K3-K5.rds") 
mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
saveRDS(mod_bf_large,file="Data/mod_bf_large_K3-K5.rds")
```

## H2: Location

```{r TestingH2,eval=TRUE,warning=FALSE,message=FALSE}
### Covid PA vs Covid US [T1-T2]
#**Difference between Mustafa T1 and Mustafa T2**
mdf12<-subset(mdf,treatment==1|treatment==2)
mdf12$treatment_name<-as.factor(mdf12$treatment_name)
mdf12$treatment_name<-relevel(mdf12$treatment_name, ref="COVID-PA")
tmp<-t.test(click~treatment_name,mdf12)
test_df[4,c("Group1","Group2")]<-levels(mdf12$treatment_name)
test_df[4,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[4]<-tmp$estimate[1]-tmp$estimate[2]
test_df[4,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[4]<-"Mustafa"
test_df$p[4]<-tmp$p.value
### Covid US vs Covid Lancaster [T1-T3]
#**Difference between Mustafa T1 and Mustafa T3**
mdf13<-subset(mdf,treatment==1|treatment==3)
mdf13$treatment_name<-as.factor(mdf13$treatment_name)
mdf13$treatment_name<-relevel(mdf13$treatment_name, ref="COVID-Lancaster")
tmp<-t.test(click~treatment_name,mdf13)
test_df[5,c("Group1","Group2")]<-levels(mdf13$treatment_name)
test_df[5,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[5]<-tmp$estimate[1]-tmp$estimate[2]
test_df[5,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[5]<-"Mustafa"
test_df$p[5]<-tmp$p.value
### Covid US vs Covid no place [T1-T4]
#**Difference between Mustafa T1 and Mustafa T4**
mdf14<-subset(mdf,treatment==1|treatment==4)
mdf14$treatment_name<-as.factor(mdf14$treatment_name)
mdf14$treatment_name<-relevel(mdf14$treatment_name, ref="COVID-US")
tmp<-t.test(click~treatment_name,mdf14)
test_df[6,c("Group1","Group2")]<-levels(mdf14$treatment_name)
test_df[6,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[6]<-tmp$estimate[1]-tmp$estimate[2]
test_df[6,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[6]<-"Mustafa"
test_df$p[6]<-tmp$p.value

### Covid PA vs Covid Lancaster [T2-T3]
#**Difference between Mustafa T2 and Mustafa T3**
mdf23<-subset(mdf,treatment==2|treatment==3)
mdf23$treatment_name<-as.factor(mdf23$treatment_name)
mdf23$treatment_name<-relevel(mdf23$treatment_name, ref="COVID-Lancaster")
tmp<-t.test(click~treatment_name,mdf23) 
test_df[7,c("Group1","Group2")]<-levels(mdf23$treatment_name)
test_df[7,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[7]<-tmp$estimate[1]-tmp$estimate[2]
test_df[7,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[7]<-"Mustafa"
test_df$p[7]<-tmp$p.value

### Covid PA vs Covid no place [T2-T4]
#**Difference between Mustafa T2 and Mustafa T4**
mdf24<-subset(mdf,treatment==2|treatment==4)
mdf24$treatment_name<-as.factor(mdf24$treatment_name)
mdf24$treatment_name<-relevel(mdf24$treatment_name, ref="COVID-PA")
tmp<-t.test(click~treatment_name,mdf24)
test_df[8,c("Group1","Group2")]<-levels(mdf24$treatment_name)
test_df[8,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[8]<-tmp$estimate[1]-tmp$estimate[2]
test_df[8,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[8]<-"Mustafa"
test_df$p[8]<-tmp$p.value

### Covid Lancaster vs Covid no place [T3-T4]
#**Difference between Mustafa T3 and Mustafa T4**
mdf34<-subset(mdf,treatment==3|treatment==4)
mdf34$treatment_name<-as.factor(mdf34$treatment_name)
mdf34$treatment_name<-relevel(mdf34$treatment_name, ref="COVID-Lancaster")
tmp<-t.test(click~treatment_name,mdf34)
test_df[9,c("Group1","Group2")]<-levels(mdf34$treatment_name)
test_df[9,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[9]<-tmp$estimate[1]-tmp$estimate[2]
test_df[9,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[9]<-"Mustafa"
test_df$p[9]<-tmp$p.value
```

```{r,eval=FALSE}
### Covid PA vs Covid US [T1-T2]
#**Difference between Mustafa T1 and Mustafa T2**
#**Bayes' factor**
#Given the observed experimental data, has the null hypothesis of no difference between Covid-PA and Covid-US or absence of an effect become more, or less credible?
#We assign informative priors assuming a linear regression model:
#$click_i \sim \text{Normal} (\alpha + \beta Treatment_i, \sigma)$
#where we take into account of the low click rates for FB ads generally as well as our assumed treatment effect sizes (small around 0.1) and consider that approximately $\beta\in(-0.5,0.5)$. We therefore set the prior mean at the midpoint of the interval and allow for some variation around either side ($\sigma=0.2$).
#test parameter with Bayes factor with simple Bayesian linear model, prior of b_treatment
loc_prior_small<- normal(location = 0, scale = 0.02)
loc_prior_med <- normal(location = 0, scale = 0.04)
loc_prior_large <- normal(location = 0, scale = 0.16)
mod_small <- stan_glm(click ~ treatment_name, data = mdf12, prior = loc_prior_small)
mod_med <- stan_glm(click ~ treatment_name, data = mdf12, prior = loc_prior_med)
mod_large <- stan_glm(click ~ treatment_name, data = mdf12, prior = loc_prior_large)
##compute bayes factor as change from prior odds to posterior odds:
mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
saveRDS(mod_bf_small,file="Data/mod_bf_small_M1-M2.rds")
mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
saveRDS(mod_bf_med,file="Data/mod_bf_med_M1-M2.rds")
mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
saveRDS(mod_bf_large,file="Data/mod_bf_large_M1-M2.rds")
### Covid US vs Covid Lancaster [T1-T3]
#**Difference between Mustafa T1 and Mustafa T3**
#**Bayes' factor**
loc_prior_small<- normal(location = 0, scale = 0.02)
loc_prior_med <- normal(location = 0, scale = 0.04)
loc_prior_large <- normal(location = 0, scale = 0.16)
mod_small <- stan_glm(click ~ treatment_name, data = mdf13, prior = loc_prior_small)
mod_med <- stan_glm(click ~ treatment_name, data = mdf13, prior = loc_prior_med)
mod_large <- stan_glm(click ~ treatment_name, data = mdf13, prior = loc_prior_large)
##compute bayes factor as change from prior odds to posterior odds:
 mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
 saveRDS(mod_bf_small,file="Data/mod_bf_small_M1-M3.rds") 
 mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
 saveRDS(mod_bf_med,file="Data/mod_bf_med_M1-M3.rds") 
 mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
 saveRDS(mod_bf_large,file="Data/mod_bf_large_M1-M3.rds")

### Covid US vs Covid no place [T1-T4]
#**Difference between Mustafa T1 and Mustafa T4**
#**Bayes' factor**
loc_prior_small<- normal(location = 0, scale = 0.02)
loc_prior_med <- normal(location = 0, scale = 0.04)
loc_prior_large <- normal(location = 0, scale = 0.16)
 mod_small <- stan_glm(click ~ treatment_name, data = mdf14, prior = loc_prior_small)
 mod_med <- stan_glm(click ~ treatment_name, data = mdf14, prior = loc_prior_med)
 mod_large <- stan_glm(click ~ treatment_name, data = mdf14, prior = loc_prior_large)
# #compute bayes factor as change from prior odds to posterior odds:
 mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
 saveRDS(mod_bf_small,file="Data/mod_bf_small_M1-M4.rds") 
 mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
 saveRDS(mod_bf_med,file="Data/mod_bf_med_M1-M4.rds") 
 mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
 saveRDS(mod_bf_large,file="Data/mod_bf_large_M1-M4.rds")

### Covid PA vs Covid Lancaster [T2-T3]
#**Difference between Mustafa T2 and Mustafa T3**
#**Bayes' factor**
loc_prior_small<- normal(location = 0, scale = 0.02)
loc_prior_med <- normal(location = 0, scale = 0.04)
loc_prior_large <- normal(location = 0, scale = 0.16)
 mod_small <- stan_glm(click ~ treatment_name, data = mdf23, prior = loc_prior_small)
 mod_med <- stan_glm(click ~ treatment_name, data = mdf23, prior = loc_prior_med)
 mod_large <- stan_glm(click ~ treatment_name, data = mdf23, prior = loc_prior_large)
# #compute bayes factor as change from prior odds to posterior odds:
 mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
 saveRDS(mod_bf_small,file="Data/mod_bf_small_M2-M3.rds") 
 mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
 saveRDS(mod_bf_med,file="Data/mod_bf_med_M2-M3.rds") 
 mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
 saveRDS(mod_bf_large,file="Data/mod_bf_large_M2-M3.rds")


### Covid PA vs Covid no place [T2-T4]
#**Difference between Mustafa T2 and Mustafa T4**
#**Bayes' factor**
loc_prior_small<- normal(location = 0, scale = 0.02)
loc_prior_med <- normal(location = 0, scale = 0.04) 
loc_prior_large <- normal(location = 0, scale = 0.16)
 mod_small <- stan_glm(click ~ treatment_name, data = mdf24, prior = loc_prior_small)
 mod_med <- stan_glm(click ~ treatment_name, data = mdf24, prior = loc_prior_med)
 mod_large <- stan_glm(click ~ treatment_name, data = mdf24, prior = loc_prior_large)
# #compute bayes factor as change from prior odds to posterior odds:
 mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
 saveRDS(mod_bf_small,file="Data/mod_bf_small_M2-M4.rds") 
 mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
 saveRDS(mod_bf_med,file="Data/mod_bf_med_M2-M4.rds") 
 mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
 saveRDS(mod_bf_large,file="Data/mod_bf_large_M2-M4.rds")


### Covid Lancaster vs Covid no place [T3-T4]
#**Difference between Mustafa T3 and Mustafa T4**
#**Bayes' factor**
loc_prior_small<- normal(location = 0, scale = 0.02)
loc_prior_med <- normal(location = 0, scale = 0.04) 
loc_prior_large <- normal(location = 0, scale = 0.16)
 mod_small <- stan_glm(click ~ treatment_name, data = mdf34, prior = loc_prior_small)
 mod_med <- stan_glm(click ~ treatment_name, data = mdf34, prior = loc_prior_med)
 mod_large <- stan_glm(click ~ treatment_name, data = mdf34, prior = loc_prior_large)
# #compute bayes factor as change from prior odds to posterior odds:
 mod_bf_small <- bayesfactor_parameters(mod_small, null = 0)
 saveRDS(mod_bf_small,file="Data/mod_bf_small_M3-M4.rds") 
 mod_bf_med <- bayesfactor_parameters(mod_med, null = 0)
 saveRDS(mod_bf_med,file="Data/mod_bf_med_M3-M4.rds") 
 mod_bf_large <- bayesfactor_parameters(mod_large, null = 0)
 saveRDS(mod_bf_large,file="Data/mod_bf_large_M3-M4.rds")
```

## H3: Profile type Refugee, Immigrant or Neither

```{r TestingH3,eval=TRUE,warning=FALSE,message=FALSE}
### Covid Refugee vs Covid Immigrant [T1-T2]
#**Difference between Kelli T1 and Kelli T2**
kdf12<-subset(kdf,treatment==1|treatment==2)
kdf12$treatment_name<-as.factor(kdf12$treatment_name)
kdf12$treatment_name<-relevel(kdf12$treatment_name, ref="COVID-Refugee")
tmp<-t.test(click~treatment_name,kdf12) 
test_df[10,c("Group1","Group2")]<-levels(kdf12$treatment_name)
test_df[10,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[10]<-tmp$estimate[1]-tmp$estimate[2]
test_df[10,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[10]<-"Kelli"
test_df$p[10]<-tmp$p.value

### Covid Immigrant vs Covid Neither [T2-T3]
#**Difference between Kelli T2 and Kelli T3**
kdf23<-subset(kdf,treatment==2|treatment==3)
kdf23$treatment_name<-as.factor(kdf23$treatment_name)
kdf23$treatment_name<-relevel(kdf23$treatment_name, ref="COVID-Immigrant")
tmp<-t.test(click~treatment_name,kdf23)
test_df[11,c("Group1","Group2")]<-levels(kdf23$treatment_name)
test_df[11,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[11]<-tmp$estimate[1]-tmp$estimate[2]
test_df[11,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[11]<-"Kelli"
test_df$p[11]<-tmp$p.value

### Covid Refugee vs Covid Neither [T1-T3]
#**Difference between Kelli T1 and Kelli T3** 
kdf13<-subset(kdf,treatment==1|treatment==3)
kdf13$treatment_name<-as.factor(kdf13$treatment_name)
kdf13$treatment_name<-relevel(kdf13$treatment_name, ref="COVID-Refugee")
tmp<-t.test(click~treatment_name,kdf13) 
test_df[12,c("Group1","Group2")]<-levels(kdf13$treatment_name)
test_df[12,c("Group1Mean","Group2Mean")]<-tmp$estimate
test_df$Mean[12]<-tmp$estimate[1]-tmp$estimate[2]
test_df[12,c("Lower","Upper")]<-tmp$conf.int
test_df$Ad[12]<-"Kelli"
test_df$p[12]<-tmp$p.value

#No Covid Refugee vs No Covid Neither
kdf45<-subset(kdf,treatment==4|treatment==5)
kdf45$treatment_name<-as.factor(kdf45$treatment_name)
kdf45$treatment_name<-relevel(kdf45$treatment_name, ref="No-COVID-Refugee")
tmp<-t.test(click~treatment_name,kdf45)
test_df2<-rbind(test_df,c("H3 Refugee effect","No-COVID-Refugee","No-COVID-Neither",tmp$estimate[1],tmp$estimate[2],tmp$estimate[1]-tmp$estimate[2],tmp$conf.int,tmp$p.value,"Kelli"))
#save
saveRDS(test_df2,file="Data/test_df2.rds")
```

## Save Testing 1, 2, 3, data frames

```{r SaveTesting,eval=TRUE,message=FALSE,warning=FALSE}
saveRDS(df,file="Data/df.rds")
saveRDS(test_df,file="Data/test_df.rds")
```

# Check BH

```{r BH}
test_df2<-readRDS("Data/test_df2.rds")
alpha =0.05
#H1
test_dfH1<-subset(test_df2,Hypothesis=="H1 Covid effect")
test_dfH1<-test_dfH1[order(test_dfH1$p),]
test_dfH1$p<-as.numeric(test_dfH1$p)
test_dfH1$bh<-p.adjust(test_dfH1$p,"BH")
#H2
test_dfH2<-subset(test_df2,Hypothesis=="H2 Location effect")
test_dfH2<-test_dfH2[order(test_dfH2$p),]
test_dfH2$p<-as.numeric(test_dfH2$p)
test_dfH2$bh<-p.adjust(test_dfH2$p,"BH")
#H3
test_dfH3<-subset(test_df2,Hypothesis=="H3 Refugee effect")
test_dfH3<-test_dfH3[order(test_dfH3$p),]
test_dfH3$p<-as.numeric(test_dfH3$p)
test_dfH3$bh<-p.adjust(test_dfH3$p,"BH")
#Mustafa
test_dfM<-subset(test_df2,Ad=="Mustafa")
test_dfM<-test_dfM[order(test_dfM$p),]
test_dfM$p<-as.numeric(test_dfM$p)
test_dfM$bh<-p.adjust(test_dfM$p,"BH")
#Kelli
test_dfK<-subset(test_df2,Ad=="Kelli")
test_dfK<-test_dfK[order(test_dfK$p),]
test_dfK$p<-as.numeric(test_dfK$p)
test_dfK$bh<-p.adjust(test_dfK$p,"BH")
```


# Figure 3: Click Rates for Ads

Replicates Figure 3 "Click Rates for Ads" in main text.

```{r Fig3Click}
# Calculates mean, sd, se and IC
mdf_sum <- mdf %>%
  group_by(treatment_name) %>%
  summarise( 
    n=n(),
    mean=mean(click),
    sd=sd(click)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))
  
mdf_sum$treatment_name<-factor(mdf_sum$treatment_name, levels=c("No-COVID-US","COVID-US","COVID-PA","COVID-Lancaster","COVID-No-Place"))

mclick<-ggplot(mdf_sum) +
  geom_bar( aes(x=treatment_name, y=mean, fill=treatment_name), stat="identity", alpha=0.75) + ggtitle("Mustafa") +
  ylab("Proportion of ad clicks") + ylim(0,0.035) + xlab("") + 
  scale_fill_viridis(discrete=TRUE,alpha=0.75,begin=0.1,end=1,option="E", name="") +
  geom_linerange( aes(x=treatment_name, y=mean, ymin=mean-ic, ymax=mean+ic), colour="gray20", alpha=0.9, size=1.3) + 
   geom_text( aes(x=treatment_name, label = round(mean,4), y = mean + 0.002), position = position_dodge(0.9), vjust = 0 ) + 
  theme_bw()+ theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none")
 #mclick
 # Calculates mean, sd, se and IC
kdf_sum <- kdf %>%
  group_by(treatment_name) %>%
  summarise( 
    n=n(),
    mean=mean(click),
    sd=sd(click)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))
  
kdf_sum$treatment_name<-factor(kdf_sum$treatment_name, levels=c("COVID-Neither","COVID-Refugee","COVID-Immigrant","No-COVID-Neither","No-COVID-Refugee"))

kclick<-ggplot(kdf_sum) +
  geom_bar( aes(x=treatment_name, y=mean, fill=treatment_name), stat="identity", alpha=0.75) + ggtitle("Kelli") +
  ylab("Proportion of ad clicks") + ylim(0,0.035) + xlab("") + 
  scale_fill_viridis(discrete=TRUE,alpha=0.75,begin=0.0,end=0.7,option="A", name="") +
  geom_linerange( aes(x=treatment_name, y=mean, ymin=mean-ic, ymax=mean+ic), colour="gray20", alpha=0.9, size=1.3) + 
   geom_text( aes(x=treatment_name, label = round(mean,4), y = mean + 0.002), position = position_dodge(0.9), vjust = 0 ) + 
  theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none")
#kclick


gA <- ggplotGrob(mclick)
gB <- ggplotGrob(kclick)
grid::grid.newpage()
grid::grid.draw(cbind(gA, gB)) #now aligned
```

# Figure 4: Ad treatment effects

Replicates Figure 4 "Ad treatment effects" in main text.

```{r Fig4Effects}
test_df<-readRDS(file="Data/test_df.rds")
test_df <- test_df %>%
  arrange(Hypothesis,Mean)
test_df$Comparison<-paste(test_df$Group1,test_df$Group2,sep=" vs ")
test_df <- test_df %>%
  mutate(Comparison=factor(Comparison,Comparison))

#colors
mycol<-c("#FDE725FF","#1F968BFF","#404788FF")

# Combine H1, H2 and H3 in one plot color based on type of H
test_df$y_colors <- ifelse(test_df$Ad=="Mustafa",mycol[3],mycol[2])

pAll<- ggplot(test_df, aes(x=Comparison, y=Mean)) +
  geom_segment(aes(x=Comparison, xend=Comparison, y=Lower, yend=Upper), color="gray", size=1.5, alpha=0.75) +
  geom_point(aes(x=Comparison, y=Mean,color=Ad),  size=3, alpha=0.85) +#size=abs(test_df$Mean)
  #Specify colours
  scale_color_manual(values = test_df$y_colors) +
  #geom_point(aes(x=5, y=0.002), colour="blue") +
  geom_hline(yintercept=0.0,color="gray") +
  theme_bw() +
  coord_flip() +
  #theme(legend.position="none"
        #,axis.text.y = element_text(color = y_colors)) +
  theme(panel.background = element_rect(fill = "white", colour = "black"),
        strip.background = element_rect(fill = "white", colour = "black"),
        legend.key = element_blank(),
        legend.title = element_blank()
        #axis.text.y = element_text(color = y_colors)
        ) +
  xlab("") +
  ylab("Difference in means estimate") +
  ggtitle("Ad treatment effects") +
  facet_grid(~Hypothesis, 
             scales = "free_x", # Let the x axis vary across facets.
             space = "free_x",  # Let the width of facets vary and force all bars to have the same width.
             switch = "x")     # Move the facet labels to the bottom.
pAll
```

# Figure 5: Positive reactions to ads

Replicates Figure 5 "Positive reactions to ads" in main text.

```{r Fig5Positive, warning=FALSE,message=FALSE}
# Calculates mean, sd, se and IC
mdf_sum <- mdf %>%
  group_by(treatment_name) %>%
  summarise( 
    n=n(),
    mean=mean(positive),
    sd=sd(positive)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))
  
mdf_sum$treatment_name<-factor(mdf_sum$treatment_name, levels=c("No-COVID-US","COVID-US","COVID-PA","COVID-Lancaster","COVID-No-Place"))

mpos<-ggplot(mdf_sum) +
  geom_bar( aes(x=treatment_name, y=mean, fill=treatment_name), stat="identity", alpha=0.75) + ggtitle("Mustafa") +
  ylab("Proportion of positive reactions") + ylim(0,0.02) + xlab("") + 
  scale_fill_viridis(discrete=TRUE,alpha=0.75,begin=0.1,end=1,option="E", name="") +
  geom_linerange( aes(x=treatment_name, y=mean, ymin=mean-ic, ymax=mean+ic), colour="gray20", alpha=0.9, size=1.3) + 
   geom_text( aes(x=treatment_name, label = round(mean,4), y = mean + 0.0012), position = position_dodge(0.9), vjust = 0 ) + 
  theme_bw()+ theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none")
#mpos
 
 # Calculates mean, sd, se and IC
kdf_sum <- kdf %>%
  group_by(treatment_name) %>%
  summarise( 
    n=n(),
    mean=mean(positive),
    sd=sd(positive)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))
  
kdf_sum$treatment_name<-factor(kdf_sum$treatment_name, levels=c("COVID-Neither","COVID-Refugee","COVID-Immigrant","No-COVID-Neither","No-COVID-Refugee"))

kpos<-ggplot(kdf_sum) +
  geom_bar( aes(x=treatment_name, y=mean, fill=treatment_name), stat="identity", alpha=0.75) + ggtitle("Kelli") +
  ylab("Proportion of positive reactions") + ylim(0,0.02) + xlab("") + 
  scale_fill_viridis(discrete=TRUE,alpha=0.75,begin=0.0,end=0.7,option="A", name="") +
  geom_linerange( aes(x=treatment_name, y=mean, ymin=mean-ic, ymax=mean+ic), colour="gray20", alpha=0.9, size=1.3) + 
   geom_text( aes(x=treatment_name, label = round(mean,4), y = mean + 0.002), position = position_dodge(0.9), vjust = 0 ) + 
  theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1),legend.position="none")
 #kpos
 
gA <- ggplotGrob(mpos)
gB <- ggplotGrob(kpos)
grid::grid.newpage()
grid::grid.draw(cbind(gA, gB)) #now aligned
```

# Appendix B: also can be run via `design.rmd`.

Replicates simulations of the experimental design — submitted for Stage 1 Report.

```{r AppendixB-1, eval=FALSE}
# For hypotheses please refer to main text.
# 
# ## Experiment 1: 5 arms
# 
# * T1 = covid - US
# * T2 = covid - PA
# * T3 = covid - Lancaster
# * T4 = covid - No location
# * T5 = no covid - US
# 
# We want to learn whether there is differential support for refugee ads on Facebook. Respondents
# are randomly assigned to receive ads with refugees with information on the above five arms. Assignment to each of the five arms is with equal probabilities, and other then mention of covid and location, ads otherwise identical. We define our outcome of interest as the difference in click rates between experimental conditions.
# 
# In settings of multiple treatment arms, we could do a number of pairwise comparisons: across treatments and each treatment against control.
# 
# #### Design Declaration A
# 
# * Model: 
# 
# We specify a population of size $N$ where a unit $i$ has a potential outcome, $Y_i(Z=0)$, when it remains untreated and $M\, (m=1,2,\cdots,M)$ potential outcomes defined according to the treatment that it receives. The effect of each treatment on the outcome of unit $i$ is equal to the difference in the potential outcome under treatment condition m and the control condition: $Y_i(Z=m)-Y_i(Z=0)$.
# 
# * Inquiry:
# 
# We are interested in all of the pairwise comparisons between arms: $E[Y(m)-Y(m')]$, for all $m\in\{1,\cdots,M\}$.
# 
# * Data strategy:
# 
# We randomly assign $k/N$ units to each of the treatment arms.
# 
# * Answer strategy:
# 
# Take every pairwise difference in means corresponding to the specific estimand.
set.seed(123)
N <- 450000 #450K
covid_effect<-0.1
us_effect<-0.05
pa_effect<-0.075
lancaster_effect<-0.08
outcome_means <- c(covid_effect+us_effect #covid-us
                   ,covid_effect+pa_effect #covid-pa
                   ,covid_effect+lancaster_effect #covid-lancaster
                   ,covid_effect #covid-nolocation
                   ,us_effect #nocovid-us
                   )
sd_i <- 0.2
outcome_sds <- c(0, 0, 0, 0, 0)

# Population
population <- declare_population(N = N, u_1 = rnorm(N, 0, outcome_sds[1L]), 
    u_2 = rnorm(N, 0, outcome_sds[2L]), u_3 = rnorm(N, 0, outcome_sds[3L]), 
    u_4 = rnorm(N, 0, outcome_sds[4L]), u_5 = rnorm(N, 0, outcome_sds[5L]),
    u = rnorm(N) * sd_i)
# Potential outcomes
potential_outcomes <- declare_potential_outcomes(formula = Y ~ (outcome_means[1] + 
    u_1) * (Z == "1") + (outcome_means[2] + u_2) * (Z == "2") + 
    (outcome_means[3] + u_3) * (Z == "3") + (outcome_means[4] + 
    u_4) * (Z == "4") + + (outcome_means[5] + 
    u_5) * (Z == "5") + u ,  conditions = c("1", "2", "3", "4", "5"), 
    assignment_variables = Z)
# Estimands
estimand <- declare_estimands(ate_Y_2_1 = mean(Y_Z_2 - Y_Z_1), ate_Y_3_1 = mean(Y_Z_3 - 
    Y_Z_1), ate_Y_4_1 = mean(Y_Z_4 - Y_Z_1), ate_Y_5_1 = mean(Y_Z_5 - Y_Z_1),
    ate_Y_3_2 = mean(Y_Z_3 - Y_Z_2), ate_Y_4_2 = mean(Y_Z_4 - Y_Z_2), 
    ate_Y_5_2 = mean(Y_Z_5 - Y_Z_2), ate_Y_4_3 = mean(Y_Z_4 - 
    Y_Z_3), ate_Y_5_3 = mean(Y_Z_5 - Y_Z_3), ate_Y_5_4 = mean(Y_Z_5 - Y_Z_4))
# Assignment
assignment <- declare_assignment(num_arms = 5, conditions = c("1", "2", "3", 
"4","5"), assignment_variable = Z)
reveal_Y <- declare_reveal(assignment_variables = Z)
# Estimator
estimator <- declare_estimator(handler = function(data) {
    estimates <- rbind.data.frame(
      ate_Y_2_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_Y_3_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_Y_4_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_Y_5_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_Y_3_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_Y_4_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_Y_5_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_Y_4_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_Y_5_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_Y_5_4 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "4", condition2 = "5")
      )
    names(estimates)[names(estimates) == "N"] <- "N_DIM"
    estimates$estimator_label <- c("DIM (Z_2 - Z_1)", "DIM (Z_3 - Z_1)", 
    "DIM (Z_4 - Z_1)", "DIM (Z_5 - Z_1)","DIM (Z_3 - Z_2)", "DIM (Z_4 - Z_2)", "DIM (Z_5 - Z_2)",
    "DIM (Z_4 - Z_3)", "DIM (Z_5 - Z_3)", "DIM (Z_5 - Z_4)")
    estimates$estimand_label <- rownames(estimates)
    estimates$estimate <- estimates$coefficients
    estimates$term <- NULL
    return(estimates)
})
multi_arm_design <- population + potential_outcomes + assignment + 
    reveal_Y + estimand + estimator

# Diagnose Experiment 1 ad click rate:
t<-Sys.time()
diagnosis <- diagnose_design(multi_arm_design)
Sys.time()-t
saveRDS(diagnosis,file="Data/diagnosis-1.rds")
```
Print table
```{r eval=TRUE}
diagnosis<-readRDS(file="Data/diagnosis-1.rds")
dat1<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","bias","rmse","power","coverage","mean_estimate","sd_estimate","mean_se","type_s_rate","mean_estimand","n_sims")]
dat2<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","se(bias)","se(rmse)","se(power)","se(coverage)","se(mean_estimate)","se(sd_estimate)","se(mean_se)","se(type_s_rate)","se(mean_estimand)","n_sims")]
dat2$estimand_label<-NA
dat2$estimator_label<-NA
tmp_n<-nrow(dat1)+nrow(dat2)
dat<-data.frame(Estimand=rep(NA,tmp_n),Estimator=rep(NA,tmp_n)
                ,Bias=rep(NA,tmp_n),RMSE=rep(NA,tmp_n)
                ,Power=rep(NA,tmp_n),Coverage=rep(NA,tmp_n)
                ,Mean_Estimate=rep(NA,tmp_n),SD_Estimate=rep(NA,tmp_n)
                ,Mean_SE=rep(NA,tmp_n) ,Type_S_Rate=rep(NA,tmp_n)
                ,Mean_Estimand=rep(NA,tmp_n),N_Sims=rep(NA,tmp_n))

j1<-j2<-1
for(i in 1:tmp_n){
  if(i%%2==0){
  dat[i,]<-dat2[j2,]
  j2<-j2+1
  }else{
  dat[i,]<-dat1[j1,]
  j1<-j1+1
  }
}
print(xtable(dat[,1:(ncol(dat)-1)],digits=2), include.rownames=FALSE)
```

```{r AppendixB-2, eval=FALSE}
#Outcome of refugee thermometer, after clicking on ad:

#Some respondents will not have thermometer ratings because of not clicking on the ads to be routed to the surveys; we can consider this as a type of attrition/missing data. This can affect both power and bias.

#As such, we set up a design that accounts for attrition.

#### Design Declaration B

# * Model: 
# 
# We specify a model with a population $N$ that has three variables affected by treatment: response variable $R_i$, outcome (here refugee thermometer rating in the survey) $Y_i$, which is correlated with response variable through parameter $\rho$. $Y_i^{obs}$ is the measured version of $Y_i$, which is only observed when $R_i=1$. For our setting, when a respondent is willing to click on the ad and answer the survey $R_i=1$.
# 
# * Inquiry:
# 
# Here we're interested in knowing the average of all respondents' differences in treatment arm potential outcomes, all of the pairwise comparisons between arms: $E[Y(m)-Y(m')]$, for all $m\in\{1,\cdots,M\}$. But we're also interested in the average treatment effect on reporting $E[R_i(m)-R_i(m')]$ as well as the pairwise comparison between treatment arms among those who report: $E[Y_i(m)-Y_i(m')|R_i=1]$.
# 
# * Data strategy:
# 
# We randomly assign $N/k = 90,000$ units to each of the treatment arms.
# 
# * Answer strategy:
# 
# $R_i$ and $Y_i^{obs}$, take every pairwise difference in means corresponding to the specific estimand.
# 
# Experiment 1: 5 arms
# 
# * T1 = covid - US
# * T2 = covid - PA
# * T3 = covid - Lancaster
# * T4 = covid - No location
# * T5 = no covid - US

#Starting parameters
N <- 450000
a_R <- 0
#Likelihood of responding to survey after exposed to treatment arm: let covid effect on going to survey be 0.3
b1_R <- 0.5 #covid - US
b2_R <- 0.6 #covid - PA
b3_R <- 0.7 #covid - Lancaster
b4_R <- 0.4 #covid - No location (-0.1 from US)
b5_R <- 0.2 #no covid - US
a_Y <- 0
#Effect on thermometer rating after exposed to treatment arm:
b1_Y <- 0.5   #covid - US
b2_Y <- 0.6   #covid - PA
b3_Y <- 0.7   #covid - Lancaster
b4_Y <- 0.4   #covid - No location (-0.1 from US)
b5_Y <- 0.2   #no covid - US
#correl
rho <- c(0.0,0.2,0.8)

#set up
t<-Sys.time()
for(i in 1:3){
  cat("Start Design:",i,"\n")
#Population
population <- declare_population(N = N, u = rnorm(N), v=rnorm(N)
              ,u1_R = rnorm(N),u2_R = rnorm(N),u3_R = rnorm(N),u4_R = rnorm(N),u5_R = rnorm(N)
              ,u1_Y = rnorm(N,mean = rho[i] * u1_R, sd = sqrt(1 - rho[i]^2)),u2_Y = rnorm(N,mean = rho[i] * u2_R, sd = sqrt(1 - rho[i]^2))
              ,u3_Y = rnorm(N,mean = rho[i] * u3_R, sd = sqrt(1 - rho[i]^2)),u4_Y = rnorm(N,mean = rho[i] * u4_R, sd = sqrt(1 - rho[i]^2))
              ,u5_Y = rnorm(N,mean = rho[i] * u5_R, sd = sqrt(1 - rho[i]^2))
              ) #one error eqn Y; one error eqn R;  errors for each condition in R; errors for each condition in Y
#Potential outcomes 
  #R
potential_outcomes_R <- declare_potential_outcomes(
  R ~ (a_R + b1_R + u1_R)* (Z == "1") + (a_R + b2_R + u2_R)* (Z == "2")  
  + (a_R + b3_R + u3_R)* (Z == "3") + (a_R + b4_R + u4_R)* (Z == "4")
  + (a_R + b5_R + u5_R)* (Z == "5") > v, conditions = c("1", "2", "3", "4", "5"), assignment_variables = Z)
  #Y
potential_outcomes_Y <- declare_potential_outcomes(
  Y ~ (a_Y + b1_Y + u1_Y)* (Z == "1") + (a_Y + b2_Y + u2_Y)* (Z == "2")
  + (a_Y + b3_Y + u3_Y)* (Z == "3") + (a_Y + b4_Y + u4_Y)* (Z == "4")
  + (a_Y + b5_Y + u5_Y)* (Z == "5") + u, conditions = c("1", "2", "3", "4", "5"), assignment_variables = Z) #>u
#Estimands: 3 types -- ATE on R, ATE on Y, ATE on Y|R
estimand <- declare_estimands(
  #ATE on R
  ate_R_2_1 = mean(R_Z_2 - R_Z_1), ate_R_3_1 = mean(R_Z_3 - R_Z_1), ate_R_4_1 = mean(R_Z_4 - R_Z_1), ate_R_5_1 = mean(R_Z_5 - R_Z_1),
  ate_R_3_2 = mean(R_Z_3 - R_Z_2), ate_R_4_2 = mean(R_Z_4 - R_Z_2), ate_R_5_2 = mean(R_Z_5 - R_Z_2),
  ate_R_4_3 = mean(R_Z_4 - R_Z_3), ate_R_5_3 = mean(R_Z_5 - R_Z_3), ate_R_5_4 = mean(R_Z_5 - R_Z_4)
  #ATE on Y
  ,ate_Y_2_1 = mean(Y_Z_2 - Y_Z_1), ate_Y_3_1 = mean(Y_Z_3 - Y_Z_1), ate_Y_4_1 = mean(Y_Z_4 - Y_Z_1), ate_Y_5_1 = mean(Y_Z_5 - Y_Z_1),
  ate_Y_3_2 = mean(Y_Z_3 - Y_Z_2), ate_Y_4_2 = mean(Y_Z_4 - Y_Z_2), ate_Y_5_2 = mean(Y_Z_5 - Y_Z_2),
  ate_Y_4_3 = mean(Y_Z_4 - Y_Z_3), ate_Y_5_3 = mean(Y_Z_5 - Y_Z_3), ate_Y_5_4 = mean(Y_Z_5 - Y_Z_4)
  #ATE on Y|R
  ,ate_YR_2_1 = mean((Y_Z_2 - Y_Z_1)[R == 1]), ate_YR_3_1 = mean((Y_Z_3 - Y_Z_1)[R == 1])
  ,ate_YR_4_1 = mean((Y_Z_4 - Y_Z_1)[R == 1]), ate_YR_5_1 = mean((Y_Z_5 - Y_Z_1)[R == 1])
  ,ate_YR_3_2 = mean((Y_Z_3 - Y_Z_2)[R == 1]), ate_YR_4_2 = mean((Y_Z_4 - Y_Z_2)[R == 1])
  ,ate_YR_5_2 = mean((Y_Z_5 - Y_Z_2)[R == 1]), ate_YR_4_3 = mean((Y_Z_4 - Y_Z_3)[R == 1])
  ,ate_YR_5_3 = mean((Y_Z_5 - Y_Z_3)[R == 1]), ate_YR_5_4 = mean((Y_Z_5 - Y_Z_4)[R == 1])
  )
#Assignment
assignment <- declare_assignment(num_arms = 5, conditions = c("1", "2", "3", "4", "5"), assignment_variable = Z)
#Reveal/Observed: ??
reveal <- declare_reveal(outcome_variables = c("R", "Y"), assignment_variables = Z)
observed <- declare_step(Y_obs = ifelse(R, Y, NA), handler = fabricate)

#Estimator
estimator <- declare_estimator(handler = function(data) {
    estimates <- rbind.data.frame(
      #ATE on R
      ate_R_2_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_R_3_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_R_4_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_R_5_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_R_3_2 = difference_in_means(formula = R ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_R_4_2 = difference_in_means(formula = R ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_R_5_2 = difference_in_means(formula = R ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_R_4_3 = difference_in_means(formula = R ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_R_5_3 = difference_in_means(formula = R ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_R_5_4 = difference_in_means(formula = R ~ Z, data = data, condition1 = "4", condition2 = "5"),
      # ATE on Y conditional on R
      ate_YR_2_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_YR_3_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_YR_4_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_YR_5_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_YR_3_2 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_YR_4_2 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_YR_5_2 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_YR_4_3 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_YR_5_3 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_YR_5_4 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "4", condition2 = "5"),
      #ATE on Y
      ate_Y_2_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_Y_3_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_Y_4_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_Y_5_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_Y_3_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_Y_4_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_Y_5_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_Y_4_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_Y_5_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_Y_5_4 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "4", condition2 = "5")
      )
    names(estimates)[names(estimates) == "N"] <- "N_DIM"
    estimates$estimator_label <- c(
      #R
      "DIM_R (Z_2 - Z_1)", "DIM_R (Z_3 - Z_1)", "DIM_R (Z_4 - Z_1)", "DIM_R (Z_5 - Z_1)","DIM_R (Z_3 - Z_2)"
      , "DIM_R (Z_4 - Z_2)", "DIM_R (Z_5 - Z_2)", "DIM_R (Z_4 - Z_3)", "DIM_R (Z_5 - Z_3)", "DIM_R (Z_5 - Z_4)"
      #Y|R
      , "DIM_Y_obs (Z_2 - Z_1)", "DIM_Y_obs (Z_3 - Z_1)", "DIM_Y_obs (Z_4 - Z_1)", "DIM_Y_obs (Z_5 - Z_1)","DIM_Y_obs (Z_3 - Z_2)"
      , "DIM_Y_obs (Z_4 - Z_2)", "DIM_Y_obs (Z_5 - Z_2)", "DIM_Y_obs (Z_4 - Z_3)", "DIM_Y_obs (Z_5 - Z_3)", "DIM_Y_obs (Z_5 - Z_4)"
      #Y
      ,"DIM_Y (Z_2 - Z_1)", "DIM_Y (Z_3 - Z_1)", "DIM_Y (Z_4 - Z_1)", "DIM_Y (Z_5 - Z_1)","DIM_Y (Z_3 - Z_2)"
      , "DIM_Y (Z_4 - Z_2)", "DIM_Y (Z_5 - Z_2)", "DIM_Y (Z_4 - Z_3)", "DIM_Y (Z_5 - Z_3)", "DIM_Y (Z_5 - Z_4)"
      )
    estimates$estimand_label <- rownames(estimates)
    estimates$estimate <- estimates$coefficients
    estimates$term <- NULL
    return(estimates)
})
multi_arm_attrition_design <- population + potential_outcomes_R + 
    potential_outcomes_Y + assignment + reveal + observed + 
    estimand + estimator
diagnoses <- diagnose_design(multi_arm_attrition_design)
saveRDS(diagnoses,paste("Data/multi_arm_attrition_design-rho",i,".rds",sep=""))
cat("Finished Design:",i," in ", Sys.time()-t,"\n")
}
Sys.time()-t
```

Print table
```{r eval=TRUE}
# Combine and print xtable
rho1<-readRDS("Data/multi_arm_attrition_design-rho1.rds")
rho2<-readRDS("Data/multi_arm_attrition_design-rho2.rds")
rho3<-readRDS("Data/multi_arm_attrition_design-rho3.rds")

dat1<-rho1$diagnosands_df
dat1$design_label<-"rho=0.0"
dat2<-rho2$diagnosands_df
dat2$design_label<-"rho=0.2"
dat3<-rho3$diagnosands_df
dat3$design_label<-"rho=0.8"

dat<-rbind(dat1,dat2,dat3)

dat1<-dat[,c("design_label","estimand_label","estimator_label","bias","rmse","power","coverage","mean_estimate","sd_estimate","mean_se","type_s_rate","mean_estimand","n_sims")]
dat2<-dat[,c("design_label","estimand_label","estimator_label","se(bias)","se(rmse)","se(power)","se(coverage)","se(mean_estimate)","se(sd_estimate)","se(mean_se)","se(type_s_rate)","se(mean_estimand)","n_sims")]


dat2$estimand_label<-NA
dat2$estimator_label<-NA

tmp_n<-nrow(dat1)+nrow(dat2)
d<-data.frame(Design=rep(NA,tmp_n),Estimand=rep(NA,tmp_n),Estimator=rep(NA,tmp_n)
                ,Bias=rep(NA,tmp_n),RMSE=rep(NA,tmp_n)
                ,Power=rep(NA,tmp_n),Coverage=rep(NA,tmp_n)
                ,Mean_Estimate=rep(NA,tmp_n),SD_Estimate=rep(NA,tmp_n)
                ,Mean_SE=rep(NA,tmp_n) ,Type_S_Rate=rep(NA,tmp_n)
                ,Mean_Estimand=rep(NA,tmp_n),N_Sims=rep(NA,tmp_n))

j1<-j2<-1
for(i in 1:tmp_n){
  if(i%%2==0){
  d[i,]<-dat2[j2,]
  j2<-j2+1
  }else{
  d[i,]<-dat1[j1,]
  j1<-j1+1
  }
}
print(xtable(d[1:60,2:(ncol(d)-1)],digits=2), include.rownames=FALSE) #rho=0.0

print(xtable(d[61:120,2:(ncol(d)-1)],digits=2), include.rownames=FALSE) #rho=0.2

print(xtable(d[121:nrow(d),2:(ncol(d)-1)],digits=2), include.rownames=FALSE) #rho=0.8
```

```{r AppendixB-3, eval=FALSE}
#### Post-experiment checks

#Given that we had to estimate the likely number of profiles on FB who would be exposed to our treatment arms, we conduct a verification of our power calculations using the actual N sizes FB ended up reaching (roughly 25K per arm):

set.seed(123)
N <- 125000 #Each arm 25K, 5 arms=125000
covid_effect<-0.1
us_effect<-0.05
pa_effect<-0.075
lancaster_effect<-0.08
outcome_means <- c(covid_effect+us_effect #covid-us
                   ,covid_effect+pa_effect #covid-pa
                   ,covid_effect+lancaster_effect #covid-lancaster
                   ,covid_effect #covid-nolocation
                   ,us_effect #nocovid-us
                   )
sd_i <- 0.2
outcome_sds <- c(0, 0, 0, 0, 0)

# Population
population <- declare_population(N = N, u_1 = rnorm(N, 0, outcome_sds[1L]), 
    u_2 = rnorm(N, 0, outcome_sds[2L]), u_3 = rnorm(N, 0, outcome_sds[3L]), 
    u_4 = rnorm(N, 0, outcome_sds[4L]), u_5 = rnorm(N, 0, outcome_sds[5L]),
    u = rnorm(N) * sd_i)
# Potential outcomes
potential_outcomes <- declare_potential_outcomes(formula = Y ~ (outcome_means[1] + 
    u_1) * (Z == "1") + (outcome_means[2] + u_2) * (Z == "2") + 
    (outcome_means[3] + u_3) * (Z == "3") + (outcome_means[4] + 
    u_4) * (Z == "4") + + (outcome_means[5] + 
    u_5) * (Z == "5") + u ,  conditions = c("1", "2", "3", "4", "5"), 
    assignment_variables = Z)
# Estimands
estimand <- declare_estimands(ate_Y_2_1 = mean(Y_Z_2 - Y_Z_1), ate_Y_3_1 = mean(Y_Z_3 - 
    Y_Z_1), ate_Y_4_1 = mean(Y_Z_4 - Y_Z_1), ate_Y_5_1 = mean(Y_Z_5 - Y_Z_1),
    ate_Y_3_2 = mean(Y_Z_3 - Y_Z_2), ate_Y_4_2 = mean(Y_Z_4 - Y_Z_2), 
    ate_Y_5_2 = mean(Y_Z_5 - Y_Z_2), ate_Y_4_3 = mean(Y_Z_4 - 
    Y_Z_3), ate_Y_5_3 = mean(Y_Z_5 - Y_Z_3), ate_Y_5_4 = mean(Y_Z_5 - Y_Z_4))
# Assignment
assignment <- declare_assignment(num_arms = 5, conditions = c("1", "2", "3", 
"4","5"), assignment_variable = Z)
reveal_Y <- declare_reveal(assignment_variables = Z)
# Estimator
estimator <- declare_estimator(handler = function(data) {
    estimates <- rbind.data.frame(
      ate_Y_2_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_Y_3_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_Y_4_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_Y_5_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_Y_3_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_Y_4_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_Y_5_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_Y_4_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_Y_5_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_Y_5_4 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "4", condition2 = "5")
      )
    names(estimates)[names(estimates) == "N"] <- "N_DIM"
    estimates$estimator_label <- c("DIM (Z_2 - Z_1)", "DIM (Z_3 - Z_1)", 
    "DIM (Z_4 - Z_1)", "DIM (Z_5 - Z_1)","DIM (Z_3 - Z_2)", "DIM (Z_4 - Z_2)", "DIM (Z_5 - Z_2)",
    "DIM (Z_4 - Z_3)", "DIM (Z_5 - Z_3)", "DIM (Z_5 - Z_4)")
    estimates$estimand_label <- rownames(estimates)
    estimates$estimate <- estimates$coefficients
    estimates$term <- NULL
    return(estimates)
})
multi_arm_design <- population + potential_outcomes + assignment + 
    reveal_Y + estimand + estimator

# Diagnose Experiment 1 ad click rate:
t<-Sys.time()
diagnosis <- diagnose_design(multi_arm_design)
Sys.time()-t
saveRDS(diagnosis,file="Data/post-diagnosis-1.rds")
```


Print table

```{r eval=TRUE}
diagnosis<-readRDS(file="Data/post-diagnosis-1.rds")
dat1<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","bias","rmse","power","coverage","mean_estimate","sd_estimate","mean_se","type_s_rate","mean_estimand","n_sims")]
dat2<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","se(bias)","se(rmse)","se(power)","se(coverage)","se(mean_estimate)","se(sd_estimate)","se(mean_se)","se(type_s_rate)","se(mean_estimand)","n_sims")]
dat2$estimand_label<-NA
dat2$estimator_label<-NA
tmp_n<-nrow(dat1)+nrow(dat2)
dat<-data.frame(Estimand=rep(NA,tmp_n),Estimator=rep(NA,tmp_n)
                ,Bias=rep(NA,tmp_n),RMSE=rep(NA,tmp_n)
                ,Power=rep(NA,tmp_n),Coverage=rep(NA,tmp_n)
                ,Mean_Estimate=rep(NA,tmp_n),SD_Estimate=rep(NA,tmp_n)
                ,Mean_SE=rep(NA,tmp_n) ,Type_S_Rate=rep(NA,tmp_n)
                ,Mean_Estimand=rep(NA,tmp_n),N_Sims=rep(NA,tmp_n))

j1<-j2<-1
for(i in 1:tmp_n){
  if(i%%2==0){
  dat[i,]<-dat2[j2,]
  j2<-j2+1
  }else{
  dat[i,]<-dat1[j1,]
  j1<-j1+1
  }
}
print(xtable(dat[,1:(ncol(dat)-1)],digits=2), include.rownames=FALSE)
```


```{r AppendixB-4, eval=FALSE}
## Experiment 2: 5 arms
# 
# * T1 = covid - Refugee
# * T2 = no covid - Refugee
# * T3 = covid - Neither
# * T4 = no covid - Neither
# * T5 = covid - Immigrant
# 
# We want to learn whether there is differential support for refugee ads on Facebook. Respondents
# are randomly assigned to receive ads with refugees with information on the above five arms. Assignment to each of the five arms is with equal probabilities, and other then mention of covid and type of individual, ads otherwise identical. We define our outcome of interest as the difference in click rates between experimental conditions.
# 
# We'll focus on pairwise comparisons across treatments (a conservative approach given our main hypotheses will be answered with comparisons of T1-T2, T2-T4, T3-T4, T1-T5, T3-T5).
# 
# #### Design Declaration A
# 
# * Model: 
# 
# We specify a population of size $N$ where a unit $i$ has a potential outcome, $Y_i(Z=0)$, when it remains untreated and $M\, (m=1,2,\cdots,M)$ potential outcomes defined according to the treatment that it receives. The effect of each treatment on the outcome of unit $i$ is equal to the difference in the potential outcome under treatment condition m and the control condition: $Y_i(Z=m)-Y_i(Z=0)$.
# 
# * Inquiry:
# 
# We are interested in all of the pairwise comparisons between arms: $E[Y(m)-Y(m')]$, for all $m\in\{1,\cdots,M\}$.
# 
# * Data strategy:
# 
# We randomly assign $k/N$ units to each of the treatment arms.
# 
# * Answer strategy:
# 
# Take every pairwise difference in means corresponding to the specific estimand.
set.seed(123)
N <- 450000 #450K
covid_effect<-0.1 #assume same covid effect as Experiment 1
refugee_effect<-0.075 #assume refugee effect is positive and larger than immigrant
immigrant_effect<-0.025 #assume immigrant effect is positive and smaller than refugee effect
outcome_means <- c(covid_effect+refugee_effect #covid - Refugee
                   ,refugee_effect #no covid - Refugee
                   ,covid_effect#covid - Neither
                   ,0 #no covid - Neither; assume no effect of ad
                   ,covid_effect+immigrant_effect#covid - Immigrant
                   )# also assumes that there are no interaction effects
sd_i <- 0.2
outcome_sds <- c(0, 0, 0, 0, 0)

# Population
population <- declare_population(N = N, u_1 = rnorm(N, 0, outcome_sds[1L]), 
    u_2 = rnorm(N, 0, outcome_sds[2L]), u_3 = rnorm(N, 0, outcome_sds[3L]), 
    u_4 = rnorm(N, 0, outcome_sds[4L]), u_5 = rnorm(N, 0, outcome_sds[5L]),
    u = rnorm(N) * sd_i)
# Potential outcomes
potential_outcomes <- declare_potential_outcomes(formula = Y ~ (outcome_means[1] + 
    u_1) * (Z == "1") + (outcome_means[2] + u_2) * (Z == "2") + 
    (outcome_means[3] + u_3) * (Z == "3") + (outcome_means[4] + 
    u_4) * (Z == "4") + + (outcome_means[5] + 
    u_5) * (Z == "5") + u ,  conditions = c("1", "2", "3", "4", "5"), 
    assignment_variables = Z)
# Estimands
estimand <- declare_estimands(ate_Y_2_1 = mean(Y_Z_2 - Y_Z_1), ate_Y_3_1 = mean(Y_Z_3 - 
    Y_Z_1), ate_Y_4_1 = mean(Y_Z_4 - Y_Z_1), ate_Y_5_1 = mean(Y_Z_5 - Y_Z_1),
    ate_Y_3_2 = mean(Y_Z_3 - Y_Z_2), ate_Y_4_2 = mean(Y_Z_4 - Y_Z_2), 
    ate_Y_5_2 = mean(Y_Z_5 - Y_Z_2), ate_Y_4_3 = mean(Y_Z_4 - 
    Y_Z_3), ate_Y_5_3 = mean(Y_Z_5 - Y_Z_3), ate_Y_5_4 = mean(Y_Z_5 - Y_Z_4))
# Assignment
assignment <- declare_assignment(num_arms = 5, conditions = c("1", "2", "3", 
"4","5"), assignment_variable = Z)
reveal_Y <- declare_reveal(assignment_variables = Z)
# Estimator
estimator <- declare_estimator(handler = function(data) {
    estimates <- rbind.data.frame(
      ate_Y_2_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_Y_3_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_Y_4_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_Y_5_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_Y_3_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_Y_4_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_Y_5_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_Y_4_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_Y_5_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_Y_5_4 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "4", condition2 = "5")
      )
    names(estimates)[names(estimates) == "N"] <- "N_DIM"
    estimates$estimator_label <- c("DIM (Z_2 - Z_1)", "DIM (Z_3 - Z_1)", 
    "DIM (Z_4 - Z_1)", "DIM (Z_5 - Z_1)","DIM (Z_3 - Z_2)", "DIM (Z_4 - Z_2)", "DIM (Z_5 - Z_2)",
    "DIM (Z_4 - Z_3)", "DIM (Z_5 - Z_3)", "DIM (Z_5 - Z_4)")
    estimates$estimand_label <- rownames(estimates)
    estimates$estimate <- estimates$coefficients
    estimates$term <- NULL
    return(estimates)
})
multi_arm_design2 <- population + potential_outcomes + assignment + 
    reveal_Y + estimand + estimator
# Diagnose Experiment 1 ad click rate:
t<-Sys.time()
diagnosis <- diagnose_design(multi_arm_design2)
Sys.time()-t
saveRDS(diagnosis,file="diagnosis-2.rds")
```

Print table

```{r eval=TRUE}
diagnosis<-readRDS(file="Data/diagnosis-2.rds")
dat1<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","bias","rmse","power","coverage","mean_estimate","sd_estimate","mean_se","type_s_rate","mean_estimand","n_sims")]
dat2<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","se(bias)","se(rmse)","se(power)","se(coverage)","se(mean_estimate)","se(sd_estimate)","se(mean_se)","se(type_s_rate)","se(mean_estimand)","n_sims")]
dat2$estimand_label<-NA
dat2$estimator_label<-NA
tmp_n<-nrow(dat1)+nrow(dat2)
dat<-data.frame(Estimand=rep(NA,tmp_n),Estimator=rep(NA,tmp_n)
                ,Bias=rep(NA,tmp_n),RMSE=rep(NA,tmp_n)
                ,Power=rep(NA,tmp_n),Coverage=rep(NA,tmp_n)
                ,Mean_Estimate=rep(NA,tmp_n),SD_Estimate=rep(NA,tmp_n)
                ,Mean_SE=rep(NA,tmp_n) ,Type_S_Rate=rep(NA,tmp_n)
                ,Mean_Estimand=rep(NA,tmp_n),N_Sims=rep(NA,tmp_n))

j1<-j2<-1
for(i in 1:tmp_n){
  if(i%%2==0){
  dat[i,]<-dat2[j2,]
  j2<-j2+1
  }else{
  dat[i,]<-dat1[j1,]
  j1<-j1+1
  }
}
print(xtable(dat[,1:(ncol(dat)-1)],digits=2), include.rownames=FALSE)
```

```{r AppendixB-5, eval=FALSE}
#### Design Declaration B

# * Model: 
# 
# We specify a model with a population $N$ that has three variables affected by treatment: response variable $R_i$, outcome (here refugee thermometer rating in the survey) $Y_i$, which is correlated with response variable through parameter $\rho$. $Y_i^{obs}$ is the measured version of $Y_i$, which is only observed when $R_i=1$. For our setting, when a respondent is willing to click on the ad and answer the survey $R_i=1$.
# 
# * Inquiry:
# 
# Here we're interested in knowing the average of all respondents' differences in treatment arm potential outcomes, all of the pairwise comparisons between arms: $E[Y(m)-Y(m')]$, for all $m\in\{1,\cdots,M\}$. But we're also interested in the average treatment effect on reporting $E[R_i(m)-R_i(m')]$ as well as the pairwise comparison between treatment arms among those who report: $E[Y_i(m)-Y_i(m')|R_i=1]$.
# 
# * Data strategy:
# 
# We randomly assign $N/k = 90,000$ units to each of the treatment arms.
# 
# * Answer strategy:
# 
# $R_i$ and $Y_i^{obs}$, take every pairwise difference in means corresponding to the specific estimand.
# 
# Experiment 2:
# 
# * T1 = covid - Refugee
# * T2 = no covid - Refugee
# * T3 = covid - Neither
# * T4 = no covid - Neither
# * T5 = covid - Immigrant

#Starting parameters
N <- 450000
a_R <- 0
#Likelihood of responding to survey after exposed to treatment arm: let covid effect on going to survey be 0.3
b1_R <- 0.1+0.075 + 0.2 #covid - refugee
b2_R <- 0.075 + 0.2  #no covid - refugee
b3_R <- 0.1 + 0.2  #covid - neither
b4_R <- 0.02 + 0.2 #no covid - neither
b5_R <- 0.1+0.025 + 0.2  #covid - immigrant
a_Y <- 0
#Effect on thermometer rating after exposed to treatment arm:
b1_Y <- 0.1+0.075   #covid - refugee
b2_Y <- 0.075   #no covid - refugee
b3_Y <- 0.1   #covid - neither
b4_Y <- 0.02  #no covid - neither
b5_Y <- 0.1+0.025   #covid - immigrant
#correl
rho <- c(0.0,0.2,0.8)

#set up
t<-Sys.time()
for(i in 1:3){
  cat("Start Design:",i,"\n")
#Population
population <- declare_population(N = N, u = rnorm(N), v=rnorm(N)
              ,u1_R = rnorm(N),u2_R = rnorm(N),u3_R = rnorm(N),u4_R = rnorm(N),u5_R = rnorm(N)
              ,u1_Y = rnorm(N,mean = rho[i] * u1_R, sd = sqrt(1 - rho[i]^2)),u2_Y = rnorm(N,mean = rho[i] * u2_R, sd = sqrt(1 - rho[i]^2))
              ,u3_Y = rnorm(N,mean = rho[i] * u3_R, sd = sqrt(1 - rho[i]^2)),u4_Y = rnorm(N,mean = rho[i] * u4_R, sd = sqrt(1 - rho[i]^2))
              ,u5_Y = rnorm(N,mean = rho[i] * u5_R, sd = sqrt(1 - rho[i]^2))
              ) #one error eqn Y; one error eqn R;  errors for each condition in R; errors for each condition in Y
#Potential outcomes 
  #R
potential_outcomes_R <- declare_potential_outcomes(
  R ~ (a_R + b1_R + u1_R)* (Z == "1") + (a_R + b2_R + u2_R)* (Z == "2")  
  + (a_R + b3_R + u3_R)* (Z == "3") + (a_R + b4_R + u4_R)* (Z == "4")
  + (a_R + b5_R + u5_R)* (Z == "5") > v, conditions = c("1", "2", "3", "4", "5"), assignment_variables = Z)
  #Y
potential_outcomes_Y <- declare_potential_outcomes(
  Y ~ (a_Y + b1_Y + u1_Y)* (Z == "1") + (a_Y + b2_Y + u2_Y)* (Z == "2")
  + (a_Y + b3_Y + u3_Y)* (Z == "3") + (a_Y + b4_Y + u4_Y)* (Z == "4")
  + (a_Y + b5_Y + u5_Y)* (Z == "5") + u, conditions = c("1", "2", "3", "4", "5"), assignment_variables = Z) #>u
#Estimands: 3 types -- ATE on R, ATE on Y, ATE on Y|R
estimand <- declare_estimands(
  #ATE on R
  ate_R_2_1 = mean(R_Z_2 - R_Z_1), ate_R_3_1 = mean(R_Z_3 - R_Z_1), ate_R_4_1 = mean(R_Z_4 - R_Z_1), ate_R_5_1 = mean(R_Z_5 - R_Z_1),
  ate_R_3_2 = mean(R_Z_3 - R_Z_2), ate_R_4_2 = mean(R_Z_4 - R_Z_2), ate_R_5_2 = mean(R_Z_5 - R_Z_2),
  ate_R_4_3 = mean(R_Z_4 - R_Z_3), ate_R_5_3 = mean(R_Z_5 - R_Z_3), ate_R_5_4 = mean(R_Z_5 - R_Z_4)
  #ATE on Y
  ,ate_Y_2_1 = mean(Y_Z_2 - Y_Z_1), ate_Y_3_1 = mean(Y_Z_3 - Y_Z_1), ate_Y_4_1 = mean(Y_Z_4 - Y_Z_1), ate_Y_5_1 = mean(Y_Z_5 - Y_Z_1),
  ate_Y_3_2 = mean(Y_Z_3 - Y_Z_2), ate_Y_4_2 = mean(Y_Z_4 - Y_Z_2), ate_Y_5_2 = mean(Y_Z_5 - Y_Z_2),
  ate_Y_4_3 = mean(Y_Z_4 - Y_Z_3), ate_Y_5_3 = mean(Y_Z_5 - Y_Z_3), ate_Y_5_4 = mean(Y_Z_5 - Y_Z_4)
  #ATE on Y|R
  ,ate_YR_2_1 = mean((Y_Z_2 - Y_Z_1)[R == 1]), ate_YR_3_1 = mean((Y_Z_3 - Y_Z_1)[R == 1])
  ,ate_YR_4_1 = mean((Y_Z_4 - Y_Z_1)[R == 1]), ate_YR_5_1 = mean((Y_Z_5 - Y_Z_1)[R == 1])
  ,ate_YR_3_2 = mean((Y_Z_3 - Y_Z_2)[R == 1]), ate_YR_4_2 = mean((Y_Z_4 - Y_Z_2)[R == 1])
  ,ate_YR_5_2 = mean((Y_Z_5 - Y_Z_2)[R == 1]), ate_YR_4_3 = mean((Y_Z_4 - Y_Z_3)[R == 1])
  ,ate_YR_5_3 = mean((Y_Z_5 - Y_Z_3)[R == 1]), ate_YR_5_4 = mean((Y_Z_5 - Y_Z_4)[R == 1])
  )
#Assignment
assignment <- declare_assignment(num_arms = 5, conditions = c("1", "2", "3", "4", "5"), assignment_variable = Z)
#Reveal/Observed: ??
reveal <- declare_reveal(outcome_variables = c("R", "Y"), assignment_variables = Z)
observed <- declare_step(Y_obs = ifelse(R, Y, NA), handler = fabricate)

#Estimator
estimator <- declare_estimator(handler = function(data) {
    estimates <- rbind.data.frame(
      #ATE on R
      ate_R_2_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_R_3_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_R_4_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_R_5_1 = difference_in_means(formula = R ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_R_3_2 = difference_in_means(formula = R ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_R_4_2 = difference_in_means(formula = R ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_R_5_2 = difference_in_means(formula = R ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_R_4_3 = difference_in_means(formula = R ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_R_5_3 = difference_in_means(formula = R ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_R_5_4 = difference_in_means(formula = R ~ Z, data = data, condition1 = "4", condition2 = "5"),
      # ATE on Y conditional on R
      ate_YR_2_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_YR_3_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_YR_4_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_YR_5_1 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_YR_3_2 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_YR_4_2 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_YR_5_2 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_YR_4_3 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_YR_5_3 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_YR_5_4 = difference_in_means(formula = Y_obs ~ Z, data = data, condition1 = "4", condition2 = "5"),
      #ATE on Y
      ate_Y_2_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_Y_3_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_Y_4_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_Y_5_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_Y_3_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_Y_4_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_Y_5_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_Y_4_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_Y_5_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_Y_5_4 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "4", condition2 = "5")
      )
    names(estimates)[names(estimates) == "N"] <- "N_DIM"
    estimates$estimator_label <- c(
      #R
      "DIM_R (Z_2 - Z_1)", "DIM_R (Z_3 - Z_1)", "DIM_R (Z_4 - Z_1)", "DIM_R (Z_5 - Z_1)","DIM_R (Z_3 - Z_2)"
      , "DIM_R (Z_4 - Z_2)", "DIM_R (Z_5 - Z_2)", "DIM_R (Z_4 - Z_3)", "DIM_R (Z_5 - Z_3)", "DIM_R (Z_5 - Z_4)"
      #Y|R
      , "DIM_Y_obs (Z_2 - Z_1)", "DIM_Y_obs (Z_3 - Z_1)", "DIM_Y_obs (Z_4 - Z_1)", "DIM_Y_obs (Z_5 - Z_1)","DIM_Y_obs (Z_3 - Z_2)"
      , "DIM_Y_obs (Z_4 - Z_2)", "DIM_Y_obs (Z_5 - Z_2)", "DIM_Y_obs (Z_4 - Z_3)", "DIM_Y_obs (Z_5 - Z_3)", "DIM_Y_obs (Z_5 - Z_4)"
      #Y
      ,"DIM_Y (Z_2 - Z_1)", "DIM_Y (Z_3 - Z_1)", "DIM_Y (Z_4 - Z_1)", "DIM_Y (Z_5 - Z_1)","DIM_Y (Z_3 - Z_2)"
      , "DIM_Y (Z_4 - Z_2)", "DIM_Y (Z_5 - Z_2)", "DIM_Y (Z_4 - Z_3)", "DIM_Y (Z_5 - Z_3)", "DIM_Y (Z_5 - Z_4)"
      )
    estimates$estimand_label <- rownames(estimates)
    estimates$estimate <- estimates$coefficients
    estimates$term <- NULL
    return(estimates)
})
multi_arm_attrition_design <- population + potential_outcomes_R + 
    potential_outcomes_Y + assignment + reveal + observed + 
    estimand + estimator
diagnoses <- diagnose_design(multi_arm_attrition_design)
saveRDS(diagnoses,paste("Data/multi_arm_attrition_design2-rho",i,".rds",sep=""))
cat("Finished Design:",i," in ", Sys.time()-t,"\n")
}
Sys.time()-t
```

Print table
```{r eval=TRUE}
# Combine and print xtable
rho1<-readRDS("Data/multi_arm_attrition_design2-rho1.rds")
rho2<-readRDS("Data/multi_arm_attrition_design2-rho2.rds")
rho3<-readRDS("Data/multi_arm_attrition_design2-rho3.rds")
dat1<-rho1$diagnosands_df
dat1$design_label<-"rho=0.0"
dat2<-rho2$diagnosands_df
dat2$design_label<-"rho=0.2"
dat3<-rho3$diagnosands_df
dat3$design_label<-"rho=0.8"

dat<-rbind(dat1,dat2,dat3)

dat1<-dat[,c("design_label","estimand_label","estimator_label","bias","rmse","power","coverage","mean_estimate","sd_estimate","mean_se","type_s_rate","mean_estimand","n_sims")]
dat2<-dat[,c("design_label","estimand_label","estimator_label","se(bias)","se(rmse)","se(power)","se(coverage)","se(mean_estimate)","se(sd_estimate)","se(mean_se)","se(type_s_rate)","se(mean_estimand)","n_sims")]


dat2$estimand_label<-NA
dat2$estimator_label<-NA

tmp_n<-nrow(dat1)+nrow(dat2)
d<-data.frame(Design=rep(NA,tmp_n),Estimand=rep(NA,tmp_n),Estimator=rep(NA,tmp_n)
                ,Bias=rep(NA,tmp_n),RMSE=rep(NA,tmp_n)
                ,Power=rep(NA,tmp_n),Coverage=rep(NA,tmp_n)
                ,Mean_Estimate=rep(NA,tmp_n),SD_Estimate=rep(NA,tmp_n)
                ,Mean_SE=rep(NA,tmp_n) ,Type_S_Rate=rep(NA,tmp_n)
                ,Mean_Estimand=rep(NA,tmp_n),N_Sims=rep(NA,tmp_n))

j1<-j2<-1
for(i in 1:tmp_n){
  if(i%%2==0){
  d[i,]<-dat2[j2,]
  j2<-j2+1
  }else{
  d[i,]<-dat1[j1,]
  j1<-j1+1
  }
}
print(xtable(d[1:60,2:(ncol(d)-1)],digits=2), include.rownames=FALSE) #rho=0.0

print(xtable(d[61:120,2:(ncol(d)-1)],digits=2), include.rownames=FALSE) #rho=0.2

print(xtable(d[121:nrow(d),2:(ncol(d)-1)],digits=2), include.rownames=FALSE) #rho=0.8
```

```{r AppendixB-6, eval=FALSE}
#### Post-experiment checks

#Given that we had to estimate the likely number of profiles on FB who would be exposed to our treatment arms, we conduct a verification of our power calculations using the actual N sizes FB ended up reaching (roughly 45K per arm):

set.seed(123)
N <- 225000 #45K per arm * 5
covid_effect<-0.1 #assume same covid effect as Experiment 1
refugee_effect<-0.075 #assume refugee effect is positive and larger than immigrant
immigrant_effect<-0.025 #assume immigrant effect is positive and smaller than refugee effect
outcome_means <- c(covid_effect+refugee_effect #covid - Refugee
                   ,refugee_effect #no covid - Refugee
                   ,covid_effect#covid - Neither
                   ,0 #no covid - Neither; assume no effect of ad
                   ,covid_effect+immigrant_effect#covid - Immigrant
                   )# also assumes that there are no interaction effects
sd_i <- 0.2
outcome_sds <- c(0, 0, 0, 0, 0)

# Population
population <- declare_population(N = N, u_1 = rnorm(N, 0, outcome_sds[1L]), 
    u_2 = rnorm(N, 0, outcome_sds[2L]), u_3 = rnorm(N, 0, outcome_sds[3L]), 
    u_4 = rnorm(N, 0, outcome_sds[4L]), u_5 = rnorm(N, 0, outcome_sds[5L]),
    u = rnorm(N) * sd_i)
# Potential outcomes
potential_outcomes <- declare_potential_outcomes(formula = Y ~ (outcome_means[1] + 
    u_1) * (Z == "1") + (outcome_means[2] + u_2) * (Z == "2") + 
    (outcome_means[3] + u_3) * (Z == "3") + (outcome_means[4] + 
    u_4) * (Z == "4") + + (outcome_means[5] + 
    u_5) * (Z == "5") + u ,  conditions = c("1", "2", "3", "4", "5"), 
    assignment_variables = Z)
# Estimands
estimand <- declare_estimands(ate_Y_2_1 = mean(Y_Z_2 - Y_Z_1), ate_Y_3_1 = mean(Y_Z_3 - 
    Y_Z_1), ate_Y_4_1 = mean(Y_Z_4 - Y_Z_1), ate_Y_5_1 = mean(Y_Z_5 - Y_Z_1),
    ate_Y_3_2 = mean(Y_Z_3 - Y_Z_2), ate_Y_4_2 = mean(Y_Z_4 - Y_Z_2), 
    ate_Y_5_2 = mean(Y_Z_5 - Y_Z_2), ate_Y_4_3 = mean(Y_Z_4 - 
    Y_Z_3), ate_Y_5_3 = mean(Y_Z_5 - Y_Z_3), ate_Y_5_4 = mean(Y_Z_5 - Y_Z_4))
# Assignment
assignment <- declare_assignment(num_arms = 5, conditions = c("1", "2", "3", 
"4","5"), assignment_variable = Z)
reveal_Y <- declare_reveal(assignment_variables = Z)
# Estimator
estimator <- declare_estimator(handler = function(data) {
    estimates <- rbind.data.frame(
      ate_Y_2_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "2"), 
      ate_Y_3_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "3"), 
      ate_Y_4_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "4"), 
      ate_Y_5_1 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "1", condition2 = "5"), 
      ate_Y_3_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "3"), 
      ate_Y_4_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "4"), 
      ate_Y_5_2 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "2", condition2 = "5"),
      ate_Y_4_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "4"),
      ate_Y_5_3 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "3", condition2 = "5"),
      ate_Y_5_4 = difference_in_means(formula = Y ~ Z, data = data, condition1 = "4", condition2 = "5")
      )
    names(estimates)[names(estimates) == "N"] <- "N_DIM"
    estimates$estimator_label <- c("DIM (Z_2 - Z_1)", "DIM (Z_3 - Z_1)", 
    "DIM (Z_4 - Z_1)", "DIM (Z_5 - Z_1)","DIM (Z_3 - Z_2)", "DIM (Z_4 - Z_2)", "DIM (Z_5 - Z_2)",
    "DIM (Z_4 - Z_3)", "DIM (Z_5 - Z_3)", "DIM (Z_5 - Z_4)")
    estimates$estimand_label <- rownames(estimates)
    estimates$estimate <- estimates$coefficients
    estimates$term <- NULL
    return(estimates)
})
multi_arm_design2 <- population + potential_outcomes + assignment + 
    reveal_Y + estimand + estimator
# Diagnose Experiment 1 ad click rate:
t<-Sys.time()
diagnosis <- diagnose_design(multi_arm_design2)
Sys.time()-t
saveRDS(diagnosis,file="Data/post-diagnosis-2.rds")
```
Print table
```{r, eval=TRUE}
diagnosis<-readRDS(file="Data/post-diagnosis-2.rds")
dat1<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","bias","rmse","power","coverage","mean_estimate","sd_estimate","mean_se","type_s_rate","mean_estimand","n_sims")]
dat2<-diagnosis$diagnosands_df[,c("estimand_label","estimator_label","se(bias)","se(rmse)","se(power)","se(coverage)","se(mean_estimate)","se(sd_estimate)","se(mean_se)","se(type_s_rate)","se(mean_estimand)","n_sims")]
dat2$estimand_label<-NA
dat2$estimator_label<-NA
tmp_n<-nrow(dat1)+nrow(dat2)
dat<-data.frame(Estimand=rep(NA,tmp_n),Estimator=rep(NA,tmp_n)
                ,Bias=rep(NA,tmp_n),RMSE=rep(NA,tmp_n)
                ,Power=rep(NA,tmp_n),Coverage=rep(NA,tmp_n)
                ,Mean_Estimate=rep(NA,tmp_n),SD_Estimate=rep(NA,tmp_n)
                ,Mean_SE=rep(NA,tmp_n) ,Type_S_Rate=rep(NA,tmp_n)
                ,Mean_Estimand=rep(NA,tmp_n),N_Sims=rep(NA,tmp_n))

j1<-j2<-1
for(i in 1:tmp_n){
  if(i%%2==0){
  dat[i,]<-dat2[j2,]
  j2<-j2+1
  }else{
  dat[i,]<-dat1[j1,]
  j1<-j1+1
  }
}
print(xtable(dat[,1:(ncol(dat)-1)],digits=2), include.rownames=FALSE)
```


# Appendix D: Post experiment survey

Replicate Table D.15.


Data pulled from `Data/Refugee_Narratives_Survey_Data.csv`. `survey` object with main variables:

* `treatment_name`: treatment arm
* `therm_refugee`: refugee thermometer 0-100
* `gender`: gender of respondent
* `educ`: educ of respondent
* `state`: state respondent lives in (see survey for number to state correspondents, alpha)
* `employment`:
  - 1: employed, not looking for work
  - 2: employed looking for work
  - 3: unemployed not looking
  - 4: unemployed looking
* `hispanic`: 1 yes, 0 no
* `race`: respondent race
* `childinhouse`: children living in household
* `party`:
  - 1: Strong Democrat
  - 2: Democrat
  - 3: Lean D
  - 4: Ind
  - 5: Lean R
  - 6: R
  - 7: Strong R
* `trumpapproval`:
  - 1: Strong approval
  - 2: Approve
  - 3: Somewhat approve
  - 4: Somewhat disapprove
  - 5: Disapprove
  - 6: Strongly disapprove
* `news_` variables (1 if uses the source, 0 otherwise)
  - `news_fb` Facebook
  - `news_social` Other social media (Twitter/Insta)
  - `news_online` Online news 
  - `news_tv` TV
  - `news_print` newspapers/journals
  - `news_radio` radio
  - `news_other` other
* `follownews`
  - 1: very closely
  - 2: somewhat 
  - 3: a little
  - 4: not at all
* `follownews_covid` : follow news post covid
  - 1: a lot more than normal
  - 2: a little more
  - 3: same
  - 4: little less
  - 5: lot less
* `RCUSA_click`: 1 if click, 0 otherwise
* `covid`:
  - 1: major threat
  - 2: minor threat
  - 3: not a threat

```{r AppendixD}
library("readxl")
#csv version
#survey<-read.csv(file="Data/Refugee_Narratives_Survey_Data.csv")
#excel version
survey<-read_excel("Data/Refugee_Narratives_Survey_Data.xlsx")
survey$therm_refugee<-survey$`On a scale from 0 to 100, where 0 equals completely unfavorable and 100 equals completely favorab...-Feelings toward refugees`
survey$treatment_name<-ifelse(survey$Treatment=="1_1","M-COVID-US"
                          ,ifelse(survey$Treatment=="1_2","M-COVID-PA"
                          ,ifelse(survey$Treatment=="1_3","M-COVID-L"
                          ,ifelse(survey$Treatment=="1_4","M-COVID-Noplace"
                          ,ifelse(survey$Treatment=="1_5","M-NoCOVID-US"
                          ,ifelse(survey$Treatment=="2_1","K-COVID-Refugee"
                          ,ifelse(survey$Treatment=="2_2","K-COVID-Immigrant"
                          ,ifelse(survey$Treatment=="2_3","K-COVID-Neither"
                          ,ifelse(survey$Treatment=="2_4","K-NoCOVID-Refugee"
                          ,ifelse(survey$Treatment=="2_5","K-NoCOVID-Neither",NA))))))))))
survey$gender<-ifelse(survey$`What is your gender?`==1,"male",ifelse(survey$`What is your gender?`==2,"female",ifelse(survey$`What is your gender?`==3,"non-binary",NA)))
survey$educ<-survey$`What is the highest level of education you have achieved?`
survey$state<-as.factor(survey$`In which US state do you currently live?`)

survey$employment<-as.factor(survey$`Are you currently employed or unemployed?`)
survey$hispanic<-ifelse(survey$`Are you of Hispanic, Latino, or Spanish origins?`==1,1
                        ,ifelse(survey$`Are you of Hispanic, Latino, or Spanish origins?`==2,0,NA))
survey$race<-ifelse(survey$`What race do you associate yourself most closely with?`==1,"White"
              ,ifelse(survey$`What race do you associate yourself most closely with?`==2,"Black"
              ,ifelse(survey$`What race do you associate yourself most closely with?`==3,"Am Ind or Al Nat"
              ,ifelse(survey$`What race do you associate yourself most closely with?`==4,"Asian"
              ,ifelse(survey$`What race do you associate yourself most closely with?`==5,"Nat Haw or Pac Isl"
              ,ifelse(survey$`What race do you associate yourself most closely with?`==6,"Other",NA))))))

survey$childinhouse<-ifelse(survey$`Do you have children living in your household?`==1,1
                     ,ifelse(survey$`Do you have children living in your household?`==2,0,NA))

survey$party<-survey$`In general, would you describe yourself as a:`
survey$trumpapproval<-survey$`Do you approve or disapprove of the way Donald Trump is handling his job as President?`
survey$religion<-ifelse(survey$`What is your present religion, if any?`==1,"Protestant"
                 ,ifelse(survey$`What is your present religion, if any?`==2,"Roman Catholic"
                 ,ifelse(survey$`What is your present religion, if any?`==3,"Mormon"
                 ,ifelse(survey$`What is your present religion, if any?`==4,"Orthodox"
                 ,ifelse(survey$`What is your present religion, if any?`==5,"Jewish"
                 ,ifelse(survey$`What is your present religion, if any?`==6,"Muslim"
                 ,ifelse(survey$`What is your present religion, if any?`==7,"Buddhist"
                 ,ifelse(survey$`What is your present religion, if any?`==8,"Hindu"
                 ,ifelse(survey$`What is your present religion, if any?`==9,"Atheist"
                 ,ifelse(survey$`What is your present religion, if any?`==10,"Agnostic"
                 ,ifelse(survey$`What is your present religion, if any?`==11,"Something else"
                 ,ifelse(survey$`What is your present religion, if any?`==12,"Nothing",NA))))))))))))

survey$news_fb<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-Online: Facebook`
survey$news_fb<-ifelse(is.na(survey$news_fb),0,1)
survey$news_social<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-Online: Other social media (e.g., Twitter, Instagram...)`
survey$news_social<-ifelse(is.na(survey$news_social),0,1)
survey$news_online<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-Online: news website or app`
survey$news_online<-ifelse(is.na(survey$news_online),0,1)
survey$news_tv<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-TV`
survey$news_tv<-ifelse(is.na(survey$news_tv),0,1)
survey$news_print<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-Print (newspapers, journals)`
survey$news_print<-ifelse(is.na(survey$news_print),0,1)
survey$news_radio<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-Radio`
survey$news_radio<-ifelse(is.na(survey$news_radio),0,1)
survey$news_other<-survey$`From which of the below sources do you receive most of your news and information? Please select a...-Other`
survey$news_other<-ifelse(is.na(survey$news_other),0,1)
survey$follownews<-survey$`How closely do you follow news and current events?`
survey$follownews_covid<-survey$`Would you say that you have been following the news more closely than normal since the emergence...`
survey$RCUSA_click<- survey$`Refugee Council USA is a non-profit, non-partisan, non-governmental organization dedicated to pro...`           
survey$covid<-survey$`Would you say that the Coronavirus (COVID-19) outbreak is a major threat, a minor threat, or not...`
#saveRDS(survey,file=Data/survey.rds")

survey_women<-as.numeric(c(round(prop.table(table(survey$gender))[1],3),0.333,NA,0.333,1,1,0.667,1,0.5,1,0))
survey_educ<-c(round(mean(survey$educ,na.rm=TRUE),3),round(tapply(survey$educ,survey$treatment_name,mean,na.rm=TRUE),3))
survey_white<-c(round(prop.table(table(survey$race))[4],3), 0.333, "",0.667,1,0,0.6,0.5,0.75,1,"")
#modal religion: tot, table(sub_survey$religion,sub_survey$treatment_name)
survey_religion<-c("Protestant","Hindu/Nothing/Catholic","","Protestant","Protestant","","Nothing","Protestant","Atheist/Muslim/Other","Protestant","")
survey_party<-c(round(mean(survey$party,na.rm=TRUE),3),round(tapply(survey$party,survey$treatment_name,mean,na.rm=TRUE),3))
tmp<-data.frame(arm=c("pooled", sort(unique(survey$treatment_name)))
           ,n=as.numeric(c(nrow(survey),c(table(survey$treatment_name))))
           ,`PropWomen` = survey_women
           , Education = survey_educ
           , `PropWhite` = survey_white
           , `ModalReligion` = survey_religion
           , `Party` = survey_party
           ,`Therm`=as.numeric(c(round(mean(survey$therm_refugee,na.rm=TRUE),3),round(tapply(survey$therm_refugee,survey$treatment_name,mean,na.rm=TRUE),3)))
           ,Trumpapprov=c(round(mean(survey$trumpapproval,na.rm=TRUE),3),round(tapply(survey$trumpapproval,survey$treatment_name,mean,na.rm=TRUE),3))
           ,RCUSA_click=c(sum(survey$RCUSA_click,na.rm=TRUE),tapply(survey$RCUSA_click,survey$treatment_name,sum,na.rm=TRUE))
           ,Covid=as.numeric(c(round(mean(survey$covid,na.rm=TRUE),3),round(tapply(survey$covid,survey$treatment_name,mean,na.rm=TRUE),3))))
rownames(tmp)<-NULL

tmp %>%
  kable(full_width = F) %>%
  #kable_paper(full_width = F) %>%
  column_spec(2, width="2cm", color = spec_color(tmp$n, begin=0.5,end=.7, na_color="gray")) %>%
  column_spec(3, width="3cm", color = spec_color(tmp$`Therm`, na_color="gray")) %>%
  column_spec(4, width="3cm") %>%
  column_spec(5, width="3cm") %>%
  column_spec(6, width="3cm") %>%
  column_spec(7, width="3cm") %>%
  column_spec(8, width="3cm") %>%
  column_spec(9, width="3cm") %>%
  column_spec(10, width="3cm") %>%
  column_spec(11, width="3cm", color = "white",
              background = spec_color(tmp$Covid, end = 0.7))
#table
is.nan.data.frame <- function(x) do.call(cbind, lapply(x, is.nan))
tmp[is.nan(tmp)] <- ""
names(tmp)<-c("Arm","N","Women","Education","White","Religion","Party","Refugee Thermometer","Trump Approval", "Click","COVID")
tmp$Arm<-c("Pooled","COVID-Immigrant","COVID-Neither","COVID-Refugee","No COVID Neither","No COVID Refugee"
           ,"COVID Lancaster", "COVID No place", "COVID PA","COVID US", "No COVID US")
kable(tmp,format="latex",booktabs = T, linesep = "", caption = "Survey response summaries. Columns `Women` and `White` present proportions of women and white respondents within the arm. Education is from 1 to 6, with 1 referring to the least amount of education achieved and 6 the most. Religion presents modal religion category in each arm. Party can take values from 1 (strong Democract) to 7 (strong Republican). The thermometer and approval variables range from 0-100. Covid is a variable that measures whether the respondent feels COVID-19 is a major threat (1), minor threat (2) or not a threat (3). Click refers to whether the respondent clicked on reading more about Refugee Council USA.")%>%
  kable_styling(latex_options=c("striped","hold_position","scale_down"))%>%
  pack_rows("Kelli ads",start_row=2,end_row=6)%>%
  pack_rows("Mustafa ads",start_row=7,end_row=11)%>%
  row_spec(0, bold=T, color="black",background="lightblue")
```

# Appendix G: Unregistered Bayes Factors

Replicates unregistered Bayes Factors Table G.16.

```{r bayesfactors}
bf <- data.frame(Group=rep(c("H1","H2"),times=c(3,6))
                 ,Name=c("COVID-US vs No COVID-US","COVID-Refugee vs No COVID-Refugee","COVID-Neither vs No COVID-Neither"
  ,"US vs No Place","PA vs No Place","PA vs US","Lancaster vs PA"
  ,"Lancaster vs No Place","Lancaster vs US")
                 ,Small=NA,Medium=NA,Large=NA)
#read in data
arm<-c("M1-M5","K1-K4","K3-K5","M1-M4","M2-M4","M1-M2","M2-M3","M3-M4","M1-M3")
for(i in 1:9){
small<-readRDS(file=paste("Data/mod_bf_small_",arm[i],".rds",sep=""))
med<-readRDS(file=paste("Data/mod_bf_med_",arm[i],".rds",sep=""))
large<-readRDS(file=paste("Data/mod_bf_large_",arm[i],".rds",sep=""))
bf$Small[i]<-small$BF[2]
bf$Medium[i]<-med$BF[2]
bf$Large[i]<-large$BF[2]
}
#output table for Appendix
knitr::kable(bf,format="latex",booktabs = TRUE
             ,align = "llccc",digits=4
             ,caption = "H1 and H2 class test Bayes factors under Gaussian priors centered at zero with small (0.02), medium (0.04) and large variances (0.16).")

```
